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ABSTRACT 


The  ability  to  predict  the  state  of  the  ground  is  essential  to  manned  and  unmanned  vehicle  mobility  and 
personnel  movement,  as  well  as  determining  sensor  performance  for  both  military  and  civilian  activities. 
As  part  of  the  Army’s  Battlespace  Terrain  Reasoning  and  Awareness  research  program,  the  1-D  dynamic 
state  of  the  ground  model  FASST  (Fast  All-season  Soil  STrength)  was  developed.  It  calculates  the  ground’s 
moisture  content,  ice  content,  temperature,  and  freeze/thaw  profiles,  as  well  as  soil  strength  and  surface 
ice  and  snow  accumulation/depletion.  The  fundamental  operations  of  FASST  are  the  calculation  of  an 
energy  and  water  budget  that  quantifies  both  the  flow  of  heat  and  moisture  within  the  soil  and  also  the 
exchange  of  heat  and  moisture  at  all  interfaces  (ground/air  or  ground/snow;  snow/air)  using  both  meteo¬ 
rological  and  terrain  data.  FASST  is  designed  to  accommodate  a  range  of  users,  from  those  who  have 
intricate  knowledge  of  their  site  to  those  who  know  only  the  site  location.  It  allows  for  22  different  terrain 
materials,  including  asphalt,  concrete,  bedrock,  permanent  snow,  and  the  USCS  soil  types.  At  a  mini¬ 
mum,  the  only  weather  information  required  is  air  temperature. 
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q,op  total  moisture  flux  at  the  surface  (m/s) 

RCI  rating  cone  index 
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Ds  snow  depth  (cm) 

f  fraction  of  snowpack  that  is  wet 

g  gravitational  constant  (98 1  cm/s2) 
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Chapter  1 

Introduction 


1.0  Background 

The  ability  to  predict  the  state  of  the  ground  is  essential  to  vehicle  mobility,  both  manned 
and  unmanned,  and  personnel  movement,  as  well  as  detennining  sensor  performance. 
Trafficability,  or  ease  of  travel,  is  dictated  by  both  soil  strength  and  surface  friction. 
Surface  friction  decreases  in  the  presence  of  an  ice  or  snow  layer  or  when  the  top  of  the 
ground  becomes  too  wet.  Soil  strength  depends  on  the  soil  type  and  the  distribution  of 
moisture  and  ice  with  depth.  For  instance,  the  presence  of  a  thawed  layer  (wet,  low 
bearing  capacity)  overlying  a  competent  layer  of  frozen  ground  has  a  negative  impact  on 
mobility  as  motion  resistance  increases  and  traction  decreases  (Sullivan  et.  al.  1997). 
Infrared  and  radar  sensor  performance  is  detennined,  in  part,  by  the  state  of  the  ground. 
Weather-impacted  state-of-the-ground  conditions  resulting  in  a  high  degree  of  clutter  can 
degrade  sensor  performance. 

1.1  Overview 

The  fundamental  operations  of  FASST  (Fast  All-seasons  Soil  STrength)  are  the 
calculation  of  an  energy  and  water  budget  that  quantifies  both  the  flow  of  heat  and 
moisture  within  the  soil  and  also  the  exchange  of  heat  and  moisture  at  all  interfaces 
(ground/air  or  ground/snow;  snow/air).  To  do  this,  FASST  uses  up  to  nine  modules 
(Section  1.2).  If  all  the  necessary  weather  parameters  (Table  1.1.1)  are  known  from  real¬ 
time  data  entry,  a  climatic  database,  or  a  mesoscale  forecast  model,  then  the  only 
modules  used  are  those  that  determine  the  state  of  the  ground.  These  modules  predict  the 
soil  temperature  and  moisture  profiles  and  ultimately  the  thickness  of  any  frozen  and/or 
thawed  soil  layers.  Otherwise,  the  solar  radiation  flux  and  the  infrared  radiation  flux  are 
calculated  using  the  two  radiation  modules  (Module  4)  and  the  snow  depth  is  calculated 
with  the  Snow  Accretion-Depletion  Module  (Module  7).  In  both  cases,  FASST  requires 
information  on  latitude,  longitude,  slope,  aspect,  elevation,  and  ground  cover.  The 
required  information  can  be  obtained  from  a  TEM  (Terrestrial  Ecosystem  Mapping) 
geographic  information  system.  FASST  is  a  one-dimensional  model  run  in  a  distributed 
mode.  That  is,  FASST  is  then  run  for  each  region  object  in  the  Area  Of  Interest  (AOI). 
Region  objects  represent  areas  with  similar  attributes  (slope,  aspect,  soil,  and  land  cover). 


Table  1.1.1  Meteorological  information  necessary  for  predicting  the  state  of  the  ground 


Radiation  (W/m2) 

Incident  and  reflected  shortwave  radiation 
(0.3-3  //);  Downwelling  and  upwelling 
longwave  radiation  (3-50//) 

Air  temperature  (°C) 

Relative  humidity  (%) 

Wind  speed  (m/s) 

Snow  depth  (cm) 

Precipitation  (mm  water//?/*) 
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The  Snow  Accretion-Depletion  Module  is  also  used  to  predict  the  amount  of  water 
available  from  snowmelt.  If  the  ground  surface  is  frozen,  the  meltwater  is  assumed  to  run 
off  and  the  moisture  content  of  the  soil  is  unchanged.  If  there  is  a  thawed  layer,  the 
moisture  content  of  the  soil  increases,  as  the  meltwater  infiltrates  the  soil.  If  the  Snow 
Accretion-Depletion  Module  is  not  being  used  to  calculate  the  amount  of  meltwater  from 
the  snowpack,  then  the  available  water  (cm)  is  assumed  to  be  0.28  times  the  change  in 
thickness  of  the  snow  cover  (cm).  A  factor  of  0.28  instead  of  0.33  is  used  to  allow  for  the 
likelihood  that  some  of  the  water  in  the  snowpack  evaporates  and  some  of  it  runs  off  once 
the  thawed  soil  layer  is  near  saturation. 


1.2  Basic  Program  Format 

The  individual  models  comprising  FASST  are  grouped  into  nine  modules.  These  modules 
accomplish  the  following  tasks: 

1.  Read  in  the  meteorological  data  and  make  any  necessary  unit  conversions  and 
assumptions.  (Chapter  9) 

2.  If  the  solar  or  infrared  flux  information  is  missing,  generate  the  appropriate  values 
based  on  the  cloud  amount  and  type  for  all  required  flux  information  except  the 
ground  emitted  IR  flux,  which  is  dependent  on  the  surface  temperature.  (Chapters 
2  and  3) 

3.  Read  in  the  control  file  containing  information  concerning  initial  conditions,  soil 
profile,  and  meteorological  data.  (Chapter  9) 

4.  Initialize  the  soil  profile  and  the  state  of  the  ground.  (Chapter  10) 

START  OF  MAIN  CALCULATION  LOOP 

5.  Calculate  the  emitted  and  net  IR  fluxes.  (Chapter  3) 

6.  Calculate  the  soil  temperature  and  volumetric  moisture  content  profiles.  (Chapters 
4  and  5) 

7.  Check  for  freezing  and/or  thawing.  Calculate  the  soil  strength.  (Chapters  5  and  6) 

8.  If  there  is  snow  or  ice  on  the  ground,  check  to  see  if  it  is  melting  or  accumulating. 
If  it  is  melting,  calculate  the  corresponding  amount  of  runoff.  (Chapters  7  and  8) 

9.  Check  the  calculation  time  increment,  update  as  necessary. 

END  OF  MAIN  CALCULATION  LOOP 

10.  Output  the  results. 
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A  graphic  representation  of  the  above  process  is  shown  in  Figure  1.1. 


Figure  1.1  Flowchart  of  the  FASST  Modules. 


1.3  References 

Sullivan,  P.M.,  C.D.  Bullock,  N.A.  Renfroe,  M.R.  Albert,  G.G.  Koenig,  L.  Peck,  and 
K.  O’Neill  (1997)  Soil  Moisture  Strength  Prediction  Model  Version  II  (SMSP  II).  U.S. 
Army  Corps  of  Engineers  Waterways  Experiment  Station,  Technical  Report  GL-97-15. 
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Chapter  2 

Shortwave  (Solar)  Radiation  Model 


2.0  Introduction 

An  important  component  of  the  surface  energy  budget  is  the  solar  net  radiation.  The  solar 
radiation  at  the  surface  consists  of  direct  and  diffuse  components.  The  relative  magnitude 
of  these  two  components  depends  in  a  large  part  on  the  optical  depth  of  the  clouds.  As  the 
cloud  optical  depth  increases,  the  direct  component  will  decrease,  while  the  diffuse 
component  may  increase  initially  if  the  optical  depth  is  small  enough.  As  the  optical 
depth  increases  beyond  some  threshold  value,  the  diffuse  component  will  decrease. 


2.1  Downwelling  Shortwave  Radiation 

Shapiro  (1987,  1982,  1972)  developed  a  simple  model  to  detennine  the  direct  and  diffuse 
shortwave  irradiance  at  the  surface  using  only  standard  surface  meteorological 
observations.  The  Shortwave  Radiation  Module  (SRM)  developed  for  FASST  uses  this 
model  to  detennine  the  downwelling  shortwave  irradiance.  The  module  also  uses  a 
simple  reflected  shortwave  irradiance  model  based  on  the  downwelling  shortwave 
irradiance  and  knowledge  of  the  surface  albedo.  Surface  albedo  is  parameterized  in  terms 
of  the  land  surface  type.  For  example,  the  albedo  for  grass  is  given  as  0.20  and  for  sand 
as  0.40.  For  each  land  type  there  will  be  an  albedo  value  based  on  that  land  type.  The 
Solar  Radiation  Module  requires  cloud  information  that  can  be  obtained  from 
observations,  climatology  or  cloud  conditions  generated  by  mesoscale  forecast  models 
such  as  the  Army’s  Boundary  layer  Forecast  Model  (BFM). 


The  basic  model  approach  involves  dividing  the  atmosphere  into  k  layers  and  assuming 

Rk+Tk+Ak=  1.  (2.1) 

R,  T,  and  A  are  the  reflectance,  transmission,  and  absorption  of  layer  k.  In  the  SRM  R,  T, 
and  A  have  been  parameterized  in  terms  of  the  solar  zenith  angle,  (p0 ,  and  the  state  of  the 
atmosphere/clouds  using  the  very  extensive  SOLMET  database  (NOAA  1979).  The 
general  form  of  the  flux  equations  is 

lksi=VA+RA  {22a) 


rk  _  rj-1  jk+ 1 

sf  “  Ik+i1s T 


(2.2b) 


Layer  k  =  0  is  the  top  of  the  atmosphere.  For  operational  use,  the  atmosphere  has  been 
divided  into  three  layers  consistent  with  the  concept  of  low,  middle,  and  high  clouds.  The 
downwelling  solar  flux  at  the  ground  (bottom  of  layer  3)  is  given  as 

Ia=TxT2TJD2  (2.3) 
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A  =  d,  (dxd2  -  A,7i2  j  -  dxR2aJ f  -  Rxa  ( T2T ) 

(2.3a) 

dJ=l-RJRj_l. 

(2.3b) 

a  is  the  ground  albedo.  Assuming  the  total  transmission  A 

can  be  specified  as  the  sum  of 

the  transmission  of  the  direct  solar  component  T  "  and  the  diffuse  solar  component  T, 

Tk  can  be  written  as 

ll 

*■ 2. 
>s‘ 

+ 

*■1. 

(2.4) 

The  direct  solar  flux  component  at  the  ground  is  given  as 

jdir  rji  dir  rj-}  dir  rji  dir  j 

1  2  3  oX 

(2.5) 

and  the  diffuse  component  as 

4s  ~t 

1 

II 

(2.6) 

2 

Ioi  (=1369.3  W/m  )  is  the  shortwave  flux  at  the  top  of  the  atmosphere.  In  order  to  solve 
the  above  equations  for  the  downwelling  direct  and  diffuse  flux  it  is  necessary  to  assume 

Tdif  =R 

1  k  'V 

(2.7) 

and  therefore  Tkir  can  be  obtained  from 

V 

l 

e-4* 

ii 

(2.8) 

Shapiro  parameterized  Tk  and  Rk  in  terms  of  the  atmospheric  and  cloud  conditions  as 
follows: 


Rk=VkPk+i}-Vk)rk 

(2.9) 

Tk  =  vpk  +  0  -  vk)h 

(2.10) 

ii 

(2.11) 

pk  is  the  cloud  reflectance  for  cloud  layer  k,  rk  the  clear  sky  reflectance  for  layer  k,  rk  the 
cloud  transmission  for  layer  k,  tk  the  clear  sky  transmission  for  layer  k,  fk  the  fractional 

cloud  amount  for  layer  k,  and  IV  is  a  cloud  weighting  factor,  p,  r,  t,  z ,  and  W  are 
parameterized  in  terms  of  the  cosine  of  the  solar  zenith  angle  using  the  SOLMET  data 
set. 

Pk  =  al  +  al  +  al  cos2  Va  +  al c°s3  Vo  (2.  12a) 

rk  =  aa°k  +  aa\  +  aa\  cos2  <po  +  aa\  cos3  (po  (2. 12b) 

4  =  K  + b\  + bl  cos2  V0  +  b\  cos3  (p0  (2.12c) 

zk  =  bb°k  +  bb\  +  bbk  cos2  cpo  +  bb\  cos3  cpo.  (2. 12d) 

The  coefficients  (aPs,  aaPs,  bPs  and  bbPs)  are  parameterized  in  terms  of  the  following 
atmospheric  and  cloud  categories:  clear,  smoke  and  haze,  thin  cirrus  and  cirrostratus, 
thick  cirrus  and  cirrostratus,  altostratus  and  altocumulus,  and  low  clouds.  The  cloud 
weighting  factor  W  is  given  as 

W  =  c0+cl  cos  cp0  +  cjk  +  cjk  cos  (p0  +  c4  cos2  (pQ  +  cJl .  (2.13) 
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The  c’s  are  parameterized  in  terms  of  the  following  cloud  categories:  thin  cirrus  and 
cirrostratus,  thick  cirrus  and  cirrostratus,  altostratus  and  altocumulus,  and  low  clouds. 
The  value  of  the  coefficients  can  be  found  in  the  code  or  in  Shapiro  (1987). 


2.2  Reflected  Shortwave  Radiation 

The  reflected  shortwave  flux  is  calculated  from  the  downwelling  shortwave  flux  and  the 
surface  albedo  as  follows 

I  A  =^si-  (2-14) 

The  surface  albedo  a  has  been  parameterized  in  terms  of  the  land  cover  type. 


2.3  Net  Shortwave  Radiation 

The  net  shortwave  radiation  on  a  horizontal  surface  is  given  as 


(2.15) 


2.4  Sloping  Surface 

All  fluxes  calculated  above  pertain  to  a  flat  surface.  If  the  surface  is  sloping,  correction 
factors  must  be  used  to  get  the  insolation  at  the  surface.  The  correction  factor  differs  for 
the  different  components  of  the  solar  flux:  downwelling  direct,  diffuse,  total,  and  the 
reflected  flux.  For  the  direct  beam,  the  corrected  downwelling  shortwave  flux  impinging 
on  a  sloping  surface  is  given  as 

I'fj  =  IdJ  [cos  (p  +  sin  (po  sin  cp  cos  «9r  /  cos  (po  \ .  (2.1 6) 

cp  is  the  slope  of  the  surface,  and  9r  is  the  relative  azimuth  and  is  defined  as  the 
difference  ( &  —  90 )  between  the  azimuth  of  the  surface  &  relative  to  north  and  the  solar 
azimuth  9o  relative  to  north.  The  corrected  diffuse  component  for  a  sloping  surface  is 
given  as 

’7  =  C  {[(! 80  -  )  (/,  +  cos  9, )  +  5,  (1  +  /,  cos  5,  )]/l  80  (1  +  /, )} .  (2.17) 

fs  is  given  as 

fs  =l  +  0.5sin^o  +2sin(2^>0).  (2.18) 

The  corrected  reflected  solar  flux  for  a  sloping  surface  is  given  as 

rt=/t(l-cosp)/2.  (2.19) 
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2.5  Net  Solar  Flux  on  a  Sloping  Surface 


The  net  solar  flux  on  a  sloping  surface  is  given  as 


(2.20) 


2.6  Shortwave  Solar  Flux  Model  Validation 

The  Smart  Weapons  Operability  Enhancement  (SWOE)  field  program  data  were  used  to 
evaluate  the  shortwave  solar  flux  model.  The  SWOE  field  programs  were  conducted 
during  the  fall  of  1992  (Grayling  1-Grayling,  MI),  the  spring/summer  of  1993  (Yuma, 
AZ),  and  the  winter/spring  of  1994  (Grayling  II-Grayling,  MI).  Each  program  ran  for 
approximately  44  days  and  meteorological  information  was  collected  at  several  locations 
in  a  2-km-square  area  with  a  temporal  frequency  of  one  minute.  Hourly  average  values  of 
the  meteorological  information  were  used  to  evaluate  the  solar  flux  model.  The  hourly 
average  values  were  computed  using  the  one-minute  values  for  a  period  of  30  minutes 
before  and  after  the  hour.  The  exception  was  the  cloud  amount  information.  During 
Grayling  I,  observations  of  the  cloud  cover  were  made  on  the  hour  from  approximately 
0700  to  1900  local.  It  should  be  realized  that  the  cloud  amount  might  not  be  consistent 
with  the  hourly  average  values  of  the  measured  solar  flux.  Since  the  exact  time  of  the 
cloud  observations  is  unknown,  it  is  not  possible  to  use  the  one-minute  observation  that 
corresponds  to  the  cloud  amount.  Cloud  amount  can  have  a  significant  impact  on  both  the 
direct  and  the  diffuse  solar  flux.  In  general,  an  increase  in  the  cloud  amount  will  decrease 
the  direct  component  of  the  solar  flux  and  will  increase  the  diffuse  (more  of  the  direct  is 
scattered  as  diffuse)  up  to  some  value  of  cloud  optical  depth.  Beyond  this  value  the 
diffuse  will  decrease.  The  analysis  presented  deals  only  with  the  daytime  solar  flux 
values,  that  is,  observations  and  model  calculated  values  when  the  solar  zenith  angle  is 
less  than  90  degrees. 

Figure  2.1  is  a  comparison  of  the  model  calculated  total  solar  flux  and  the  corresponding 
measured  flux  for  a  10-day  period.  Some  of  the  differences  can  be  directly  attributed  to 
the  fact  that  the  measured  values  are  one-hour  averages  while  the  calculated  values 
basically  use  the  instantaneous  cloud  observation.  The  phasing  and  the  amplitudes  match 
fairly  well  for  clear  conditions.  Differences  occur  mainly  for  cloudy  conditions.  The 
components  of  the  solar  flux  are  presented  in  Figures  2.2  and  2.3.  As  the  cloud  amount 
increases,  the  direct  solar  flux  decreases  while  the  diffuse  tends  to  increase.  The  solar 
flux  model  does  not  require  cloud  optical  depth  information.  The  dependence  on  optical 
depth  is  implicit  in  the  coefficients  for  transmission,  reflectance,  and  the  cloud  weighting 
function.  As  indicated  above,  the  diffuse  component  should  decrease  when  the  optical 
depth  exceeds  some  critical  value. 
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Time  Series  of  the  Total  Measured  &  Calculated  Solar  Flux 


Figure  2.1  Comparison  of  the  measured  and  calculated  total  solar  flux  for  a  10-day 
period  for  Grayling  I.  Only  the  values  of  the  solar  flux  when  the  zenith  angle  is  less  than 
90  degrees  are  plotted. 


Time  Series  of  the  Direct  Measured  &  Calculated  Solar  Flux 


Figure  2.2  Same  as  Figure  2. 1  but  for  the  direct  component  of  the  solar  flux. 


8 


Time  Series  of  the  Diffuse  Measured  &  Calculated  Solar  Flux 


Figure  2.3  Same  as  Figure  2.1  but  for  the  diffuse  solar  flux. 


Most  surface  energy  budget  models,  including  FASST,  utilize  only  the  total  solar  flux. 
Figure  2.4  is  a  scatter  diagram  of  the  total  measured  and  calculated  solar  flux.  A  least 
square  fit  to  the  data  and  its  correlation  coefficient,  R2,  is  presented  in  the  figure.  The 
least  square  fit  has  been  forced  through  the  origin  (0,0). 


Grayling  I  Total  Solar  Flux 


Figure  2.4  Scatter  diagram  of  the  total  solar  measured  and  calculated  flux  and  the  least 
square  fit  to  the  data. 
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The  amount  of  scatter,  relative  to  the  least  square  fit,  increases  in  the  flux  range  from 
approximately  50  to  400  W/m  .  Some  of  this  scatter  is  no  doubt  due  to  the  comparison  of 
the  average  flux  conditions  with  calculated  fluxes  based  on  the  instantaneous  observation 
of  cloud  amount  on  the  hour.  Table  2.6.1  presents  the  equation  for  the  least  square  fit  and 
its  correlation  coefficient  for  the  total  solar  flux  and  each  of  the  flux  components. 


Table  2.6.1  Equation  of  fit  (forced  through  the  origin)  and  correlation  coefficient  for  Grayling  I. 


Solar  Flux 

Equation  of  Fit 

R2 

Total 

y=  1.036  lx 

0.8049 

Direct 

v=1.1353x 

0.8631 

Diffuse 

y= 0.8757x 

0.3752 

The  comparison  of  the  diffuse  components  shows  considerable  amount  of  scatter  as 
evident  by  the  low  correlation  coefficient.  There  are  a  number  of  data  points  where  the 
calculated  is  on  the  order  of  100  to  200  W/m2  greater  than  the  measured  value.  This 
relatively  large  difference  is  due  in  part  to  inappropriate  cloud  representation.  This 
includes  both  the  cloud  amounts  and  the  implicit  cloud  optical  depth  as  discussed  earlier. 

Figure  2.5  explores  the  impact  of  cloud  conditions  on  the  difference  of  the  calculated 
total  solar  flux  minus  the  total  solar  measured  flux.  Some  interesting  features  are 
observed  in  Figure  2.5.  The  range  of  flux  differences  is  smallest  when  there  is  no  low 
cloud  (clear,  middle,  high,  or  middle  plus  high).  When  there  are  only  low  clouds  the 
differences  are  evenly  distributed  from  -300  to  +300  W/m 2 .  When  there  are  low  clouds 
plus  another  cloud  type  (low  +  middle,  low  +  high,  or  low  +  middle  +  high),  the 
calculated  flux  tends  to  be  greater  than  the  measured  flux.  The  cloud  optical  depth  is  a 
function  of  the  cloud  thickness.  The  cloud  model  does  not  directly  use  cloud  optical 
depth  or  cloud  thickness.  One  of  the  reasons  for  this  is  that  cloud  thickness  and  cloud 
optical  depth  are  not  reported  meteorological  variables. 
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Calculated  Minus  the  Measured  Total  Solar  Flux 
for  Different  Cloud  Conditions 


Figure  2.5  Calculated  minus  the  measured  total  solar  flux  for  different  cloud  conditions 
for  Grayling  I.  The  dashed  lines  are  ±  50  W/m“. 

The  analysis  of  the  Yuma  data  indicates  less  scatter  (Figure  2.6).  This  is  anticipated, 
since  there  are  fewer  observations  with  clouds.  Again,  the  correlation  coefficient  (Table 
2.6.2)  is  lowest  for  the  diffuse  component  of  the  solar  flux.  In  general,  the  model 
performs  better  in  the  low  cloud  amount  environment  associated  with  the  meteorological 
conditions  at  Yuma.  The  flux  difference  as  a  function  of  cloud  conditions  (Figure  2.7) 
again  indicates  that  the  model  fluxes  tend  to  be  greater  than  the  measured  fluxes,  even  for 
clear  conditions.  Yuma  is  a  very  dry  environment  and  a  dust  layer  frequently  developed 
by  the  end  of  the  day  as  a  result  of  the  military  traffic  in  the  local  training  area.  The 
model  does  not  have  a  provision  for  handling  the  reduction  of  solar  flux  at  the  surface 
due  to  a  dust  layer. 


Table  2.6.2  Equation  of  fit  (forced  through  the  origin)  and  correlation  coefficient  for  Yuma. 


Solar  Flux 

Equation  of  Fit 

R2 

Total 

y=  1.0656x 

0.9515 

Direct 

y=1.0706x 

0.9023 

Diffuse 

y=0.9708x 

0.6191 

The  analysis  for  Grayling  II  is  presented  in  Figures  2.8  and  2.9  and  Table  2.6.3.  As  with 
Grayling  I,  the  correlation  coefficients  for  the  total,  direct,  and  diffuse  solar  flux  are  less 
than  the  values  for  the  Yuma  analysis. 


Table  2.6.3  Equation  of  fit  (forced  through  the  origin)  and  correlation  coefficient  for  Grayling  II. 


Solar  Flux 

Equation  of  Fit 

R 2 

Total 

v=1.0395x 

0.8363 

Direct 

y=0.98x 

0.8266 

Diffuse 

y=1.1477x 

0.5471 
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Yuma  Total  Solar 


Figure  2.6  Scatter  diagram  of  the  total  solar  measured  and  calculated  flux  and  the  least 
square  fit  to  the  data. 


Calculated  Minus  the  Measured  Total  Solar  Flux 
for  Different  Cloud  Conditions 


Figure  2.7  Calculated  minus  the  measured  total  solar  flux  for  different  cloud  conditions 
for  Yuma.  The  dashed  lines  are  ±  50  W/ni . 
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Grayling  II  Total  Solar  Flux 


Figure  2.8  Scatter  diagram  of  the  total  solar  measured  and  calculated  flux  and  the  least 
square  fit  to  the  data. 


600  - 

500  - 

400  - 

300  - 


Calculated  Minus  the  Measured  Total  Solar  Flux 
for  Different  Cloud  Conditions 


Figure  2.9  Calculated  minus  the  measured  total  solar  flux  for  different  cloud  conditions 
for  Grayling  II.  The  dashed  lines  are  ±  50  W/nf . 
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2.7  Module  Input  and  Output 


The  solar  model  inputs  are 

1 .  The  latitude  and  longitude  of  the  location. 

2.  The  Day-Of-the-Year  (DOY)  and  hour  of  the  day. 

3.  The  land  surface  category. 

4.  Low,  middle,  and  high  fractional  cloud  amount.  If  all  three  cloud  amounts  are 
missing,  the  low  cloud  amount  is  set  equal  to  the  global  average  cloud  amount  (0.54). 

5.  The  low  middle,  and  high  cloud  type.  If  missing,  the  model  uses  a  set  of  default 
conditions.  The  low  cloud  default  is  cumulus,  the  middle  cloud  type  is  altostratus,  and 
the  high  type  is  cirrus. 

Items  1  and  2  are  obtained  from  the  meteorological  data  file,  as  are  items  4  and  5  if  they 
are  available.  Item  3  is  found  in  the  FASST  input  file. 

The  model  outputs  are  the  downwelling,  reflected  and  net  shortwave  flux  (W/m  )  at  the 
surface  for  the  input  conditions.  If  the  input  conditions  are  representative  of  some  time 
period,  then  the  number  of  Joules/m  for  that  period  is  given  as  the  flux  ( W/m )  times  the 
length  of  the  period  in  seconds. 
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Chapter  3 

Longwave  (IR)  Radiation  Model 


3.0  Introduction 

The  net  longwave  radiation  at  the  surface  is  an  important  component  of  any  surface 
energy  budget  model.  Unfortunately,  measurements  of  the  net  longwave  radiation 
(upwelling  longwave  radiation  and  downwelling  longwave  radiation)  are  not  readily 
available.  Therefore,  it  is  necessary  to  use  parameters  that  are  available  and  simple 
models  to  calculate  the  upwelling  and  downwelling  longwave  radiation  at  the  surface. 
The  net  longwave  radiation  is  a  function  of  atmospheric  conditions,  especially  cloud 
conditions,  and  the  surface  optical  and  physical  properties. 


3.1  Downwelling  Longwave  Radiation 

One  of  the  major  components  to  the  downwelling  radiation  at  the  surface  is  the 
atmosphere.  Most  of  the  downwelling  atmospheric  radiation  originates  in  the  lower 
atmosphere  (below  altitudes  of  several  kilometers)  where  the  absolute  humidity  is 
relatively  high.  The  component  ( )  for  a  clear  sky  is  calculated  from 

ltl=WTa  (3-1) 

where  sa  is  the  effective  atmospheric  emissivity,  Ta  (K)  is  the  ambient  air  temperature, 

and  cr  (5.669  x  10  8  W/m2 ■  K4)  is  the  Stefan-Boltzmann  constant.  The  effective 
atmospheric  emissivity  is  calculated  from  Crawford  et  al.  (1999) 

sa  =  \24*(eJT,)m.  (3.2) 

ea  is  the  vapor  pressure  in  millibars  and  can  be  found  from 


RH 


100-^ 


(3.3) 


where  RH  is  the  relative  humidity  in  percent  and 


C,  =  Co  exP 


"  / 

(  l 

i  y 

K 

Jo  ' 

~Ta)_ 

(3.4) 


is  the  saturated  vapor  pressure.  eso  is  the  saturated  vapor  pressure  at  standard  conditions 
(«6.1 1  mbars),  l  is  the  temperature-dependent  latent  heat  of  condensation  (sublimation), 
Rv  is  the  gas  constant  for  water  vapor,  and  7V  is  the  standard  temperature  (273.15  °C). 


For  cloudy  skies  the  downwelling  flux  ( /"['  )  is  given  as 

i;!i=cIxl+c:ffxm+cheffXh 


(3.5) 
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where  c,  is  the  fractional  low  cloud  amount,  c'f  and  cheff  are  the  effective  middle  and  high 
cloud  cover  amounts  based  on  the  principle  of  random  overlap,  and  Zu,m,h)  's  the  cloud 
irradiance  factor  for  low,  middle,  and  high  clouds,  respectively.  c"'f  and  cheff  are  given  as 

<%r  =  cm(  1-c,)  and  chejf  =  (3.6) 

where  cm  (c/,)  are  the  fractional  middle  (high)  cloud  cover  amounts.  Hodges  et  al.  (1983) 
originally  defined  zu,m,h)  as 

X(i,m,h)  =  80  —  5  Z{lmh^  (3.7) 

where  Z{l  m  h)  is  either  the  low,  middle,  or  high  cloud  base  altitude  in  kilometers.  Equation 

(3.7)  has  been  modified  based  on  the  Geophysics  Laboratory  model  atmospheres  (tropics, 
mid  latitude  summer,  mid-latitude  winter,  subarctic  summer,  and  subarctic  winter)  and  is 
given  as 

%u,m,h )  =  94  —  5.8  Z{lmh) .  (3.8) 

The  cloud  base  altitude,  if  not  available  from  observations,  has  been  parameterized  in 
tenns  of  season  and  latitude  following  the  approach  by  Stowe  et  al.  (1980)  and  London 
(1957)  and  data  from  the  Global  Distribution  of  Total  Cloud  Cover  and  Cloud  Type 
Amounts  Over  Land  (1986)  prepared  by  DOE  and  NCAR  (Table  3.1.0).  The  basic 
equation  is  given  as 

^ \l,m,h )  =a(l,m,h)  (3-9) 

where  X  is  the  latitude  and  a,im,h),  b(iy/nj,),  C(ijn,h)  and  d(i,m,h)  are  parameterization 
coefficients  given  in  Table  3.1.0. 


Table  3.1.0  Coefficients  for  determining  cloud  base  altitude  ( meters )  as  a  function  of  season, 
latitude,  and  cloud  type. _ 


Winter  (Dec,  Jan,  and  Feb) 

Cloud  Type 

Latitude 

a 

b 

c 

d 

low 

>=25 

1050 

600 

1.5 

25 

low 

<25 

1050 

600 

5.0 

25 

middle 

>=25 

4100 

2000 

1.7 

25 

<25 

4100 

300 

4 

25 

high 

all 

7000 

1500 

3 

30 

Non-Winter  (Mar,  Apr,  May,  Jun,  Jul,  Aug,  Sej 

p,  Oct,  and  Nov) 

Cloud  Type 

Latitude 

a 

b 

c 

d 

low 

>=25 

1150 

600 

1.5 

25 

low 

<25 

1150 

450 

5.0 

25 

middle 

>=25 

4400 

1200 

3.0 

25 

middle 

<25 

4400 

300 

4 

25 

high 

all 

7000 

1500 

3 

30 

The  total  downwelling  radiation  is  given  as 

j  jclr  .  jcld 

ir-l  irl  irl 


(3.10) 
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3.1.1  Modification  of  the  Cloud  Irradiance  Factor  (CIF) 


The  total  downwelling  radiation  at  the  ground,  assuming  a  single  optically  thick  overcast 
cloud  layer,  can  be  written  as 

I,rl=e,aT:+(iO-5Z)  (3.11) 

based  on  the  original  CIF  (Z  is  the  cloud  base  height).  The  total  downwelling  radiation 
can  also  be  written  as 

Iirl=Sa(jTa  +(y-£a)°Tc  ■  (3-12) 


Cloud  Irradiance 


Mid-Latitude  Summer  Atmosphere 


^Calculated  -*■  Least  Square  Fit 


Figure  3.1  Calculated  cloud  irradiance  in  accordance  with  Equation  (3.12)  and  the  least 
square  fit  to  the  data  points.  Calculations  are  based  on  a  mid-latitude  summer 
atmosphere. 

The  term  in  parentheses  is  the  atmospheric  transmission  below  the  cloud  and  Tc  (K)  is 

the  cloud  base  temperature.  The  cloud  emissivity  is  assumed  to  be  one,  implying  the 
cloud  reflection  and  transmission  are  zero.  Using  a  cloud  emissivity  of  one  and  the  cloud 
base  temperature  (clouds  radiate  at  a  temperature  slightly  colder  than  the  cloud  base 
temperature)  will  result  in  an  overestimation  of  the  emitted  cloud  radiance.  This  is  offset 
by  the  fact  that  the  cloud  reflection  is  assumed  to  be  zero  and  therefore  the  cloud  does  not 
reflect  downward  any  of  the  earth-emitted  radiance.  Equating  Equations  (3.11)  and 
(3.12),  we  can  re-derive  the  term  in  parentheses  in  Equation  (3.11)  by  computing  the 
cloud  contribution  to  the  total  downwelling  radiation  using  either  measured  moisture  and 
atmospheric  temperature  profiles  or  the  Geophysics  Laboratory  model  atmospheres.  The 
atmospheric  emissivity  was  calculated  using  the  Geophysics  Laboratory  model 
atmospheres  and  Equation  (3.2)  (and  Equation  [3.3]  etc.).  Next,  the  cloud  radiance  was 
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calculated  for  overcast  optically  thick  clouds  at  altitudes  from  1  km  to  10  km  for  each  of 
the  Geophysics  Laboratory  model  atmospheres.  An  example  of  the  resulting  calculations 
is  given  in  Figure  3.1. 

The  results  for  all  Geophysics  Laboratory  model  atmospheres  and  the  correlation 
coefficients  are  given  in  Table  3.1.1. 


Table  3.1.1  Calculated  regressions  for  determining  the  cloud  irradiance  for  the  Geophysics 
Laboratory  Model  Atmospheres. _ _ 


Atmosphere 

Cloud  Irradiance 

Correlation  coefficient 

Tropical 

95-5. 8Z 

0.996 

Mid-latitude  Summer 

96-5. 8Z 

0.997 

Mid-latitude  Winter 

93-5. 8Z 

0.997 

Subarctic  Summer 

95-5.9Z 

0.998 

Subarctic  Winter 

90-5. 8Z 

0.991 

3.2  Downwelling  Infrared  Flux  Model  Validation 

The  Smart  Weapons  Operability  Enhancement  (SWOE)  field  program  data  were  used  to 
evaluate  the  downwelling  infrared  flux  model.  The  SWOE  field  programs  were 
conducted  during  the  fall  of  1992  (Grayling  I-Grayling,  MI),  the  winter/spring  of  1994 
(Grayling  II-Grayling,  MI),  and  the  spring/summer  of  1993  (Yuma,  AZ).  Each  program 
ran  for  approximately  44  days  and  meteorological  information  was  collected  at  several 
locations  in  a  2 -/cm -square  area  with  a  temporal  frequency  of  one  minute.  Hourly  average 
values  of  the  meteorological  information  were  used  to  evaluate  the  downwelling  infrared 
model.  The  hourly  average  values  were  computed  using  the  one-minute  values  for  a 
period  of  30  minutes  before  and  after  the  hour.  The  exception  was  the  cloud  amount 
information.  During  Grayling  I,  observations  of  the  cloud  cover  were  made  on  the  hour 
from  approximately  0700  to  1900  local.  During  the  night,  the  cloud  amounts  were 
derived  from  ceilometer  measurements.  The  hourly  fractional  cloud  amount  was 
computed  from  the  number  of  one-minute  ceilometer  measurements  that  indicated 
clouds.  The  sixty  observations  were  centered  on  the  hour.  Unfortunately,  the  ceilometer 
used  to  compute  the  nighttime  cloud  amounts  was  an  older  model  that  had  an  effective 
range  of  only  several  kilometers.  It  is  highly  unlikely  that  middle  and  high  clouds  ever 
entered  into  the  calculation  of  the  cloud  amount.  Figure  3.2  is  a  plot  of  the  measured  and 
computed  downwelling  longwave  flux  as  a  function  of  the  fractional  day  of  the  year. 
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Grayling  I  Downwelling  IR  Flux  Comparison 
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Figure  3.2  Comparison  of  the  downwelling  measured  and  computed  longwave  flux  based 
on  the  Grayling  I  meteorological  observations. 

The  comparison  can  also  be  presented  as  a  scatter  diagram  (Figure  3.3).  Also  is  the  linear 
least  square  fit  and  the  R2  coefficient. 


Graylin  I  Downwelling  IR  Flux  Comparison 


Figure  3.3  Comparison  of  the  measured  and  model  calculated  downwelling  longwave 
flux. 
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Some  scatter  is  evident  in  Figure  3.3.  All  points  where  the  difference  between  the 
measured  and  calculated  downwelling  fluxes  was  greater  than  50  W/m~  were  investigated. 
Fifty-two  points  out  of  932  points  fell  in  this  category.  In  every  case,  the  measured  flux 
was  greater  than  the  calculated.  In  13  of  these  cases,  the  observer  reported  frost  on  the 
dome  of  the  instrument  during  the  morning  inspection  of  the  meteorological  instruments. 
Of  the  39  remaining  cases,  24  occurred  during  periods  of  reduced  visibility  (fog  or  haze) 
and  clear  or  scattered-cloud  sky  conditions.  During  periods  of  reduced  visibility  the 
sensor  does  not  “see”  the  cold  sky,  but  some  temperature  associated  with  the  low 
atmosphere.  The  infrared  model  will  not  handle  periods  of  reduced  visibility  effectively. 
Only  two  cases  occurred  during  the  day  with  no  obscuration  to  the  visibility  and  when  an 
observer  was  on  duty  to  take  cloud  observations.  The  average  flux  values  for  clear, 
scattered  (0.1  to  0.5  cloud  cover),  broken  (0.6  to  0.9  cloud  cover),  and  overcast  (1.0  cloud 
cover)  are  given  in  Table  3.2.1. 


Table  3.2.1  Grayling  I  average  downwelling  IR  flux  for  different  cloud  conditions. 


Cloud  amount 

Measured  IR 

Calculated  IR 

clear 

284.6 

257.3 

Scattered  (0.1  to  0.5  cloud  cover) 

293.3 

282.3 

Broken  (0.5  to  0.9  cloud  cover) 

331.9 

336.3 

Overcast  ( 1 .0  Cloud  cover) 

349.8 

359.2 

A  similar  analysis  of  the  Yuma  data  is  presented  in  Figures  3.4  and  3.5. 


Yuma  Downwelling  Infrared  Flux  Comparison 


Figure  3.4  Comparison  of  the  downwelling  measured  and  computed  longwave  flux  based 
on  the  Yuma  meteorological  observations. 
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Yuma  Downwelling  Infrared  Flux  Comparison 
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Figure  3.5  Comparison  of  the  measured  and  model  calculated  downwelling  longwave 
flux. 

There  is  considerable  spread  in  the  data  resulting  in  a  relatively  low  R2.  During  the  field 
program,  daytime  temperatures  were  as  high  as  3 1  °C  with  relative  humidity  as  low  as  3 
percent.  Out  of  1124  observations,  a  difference  of  50  W/m2  or  greater  between  the 
measured  value  and  the  corresponding  calculated  occurred  only  seventeen  times.  For  four 
of  the  seventeen  cases  the  calculated  flux  was  greater  than  the  measured  flux.  In  every 
one  of  these  cases,  eight  tenths  or  more  of  thin  cirrus  was  reported.  Two  cirrus  regimes 
invaded  this  area:  mid-latitude  cirrus  and  cirrus  associated  with  a  subtropical  jet  stream. 
Cirrus  associated  with  the  subtropical  jet  usually  consists  of  smaller  ice  crystals.  In 
addition,  the  base  of  the  cirrus  normally  occurs  at  a  higher  altitude  than  mid-latitude 
cirrus  (based  on  LIDAR  and  aircraft  observations  taken  during  the  First  ISCCP  Regional 
Experiments  [FIRE]  program  at  Coffeeville,  KS).  Thin  cirrus  (meteorologically,  thin 
cirrus  is  defined  as  cirrus  that  is  transparent  enough  that  blue  sky  can  be  observed)  will 
have  an  emissivity  that  is  less  than  optically  thick  cirrus.  Paltridge  and  Platt  (1981) 
reported  cirrus  emittance  as  a  function  of  Ice  Water  Path  (IWP).  Emittance  values  range 
from  0.1  for  IWPs  of  approximately  10  g/m2  to  0.9  for  IWPs  of  60  g/m2.  The  model 
implicitly  assumes  an  emissivity  of  one  for  all  clouds.  The  calculated  base  of  the  clouds 
was  on  the  order  of  seven  kilometers.  Using  an  emissivity  (0.5  rather  than  1.0)  and  cloud 
base  altitude  (10  km  rather  than  7  km)  that  is  more  characteristic  of  thin  subtropical  jet 
stream  cirrus,  the  contribution  to  the  flux  at  the  surface  is  reduced  by  approximately  36 
W/m2.  Six  of  the  observations  that  exceeded  the  50  W/m2  criteria  occurred  when  the 
relative  humidity  was  5  percent  or  less  and  the  temperature  was  approximately  30  °C.  It  is 


y=  0.976X  +  4.2104 
R2  =  0.6221 
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not  known  whether  the  parameterization  of  the  atmospheric  emissivity  can  effectively 
handle  these  extreme  conditions.  In  addition,  the  computed  downwelling  infrared  flux  is 
fairly  sensitive  to  the  low  cloud  amount  as  anticipated.  For  example,  a  change  of  0.1  in 
the  cloud  amount  for  a  cloud  at  one  kilometer  will  result  in  a  change  of  approximately  90 
W/rrf  in  the  downwelling  flux.  The  average  flux  values  for  clear,  scattered  (0.1  to  0.5 
cloud  cover),  broken  (0.6  to  0.9  cloud  cover),  and  overcast  (1.0  cloud  cover)  are  given  in 
Table  3.2.2. 


Table  3.2.2  Yuma  average  downwelling  IR  flux  for  different  cloud  conditions. 


Cloud  amount 

Measured  IR 

Calculated  IR 

Clear 

311.7 

296.40 

Scattered  (0.1  to  0.5  cloud  cover) 

328.7 

332.33 

Broken  (0.5  to  0.9  cloud  cover) 

369.2 

373.10 

Overcast  ( 1 .0  Cloud  cover) 

348.7 

370.40 

Grayling  II  was  a  winter  field  program.  Fog,  including  ice  fog,  and  very  low  surface 
temperatures  (-20  °C)  will  tax  the  ability  of  the  model  to  accurately  predict  the 
downwelling  infrared  flux.  The  comparison  of  the  measured  and  the  calculated  flux  is 
presented  in  the  following  figures  and  table. 


Grayling  II  Downwelling  Infrared  Flux  Comparison 


Figure  3.6  Comparison  of  the  downwelling  measured  and  computed  longwave  flux  based 
on  the  Grayling  II  meteorological  observations. 
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Grayling  II  Downwelling  Infrared  Flux  Comparison 
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Figure  3.7  Comparison  of  the  measured  and  model  calculated  downwelling  longwave 
flux. 

2 

Four  percent  (49  out  of  1024)  of  the  calculated  flux  values  differed  by  more  than  50  W/m 
from  the  corresponding  measured  values.  In  25  of  the  49  cases,  either  ice/snow  was 
reported  on  the  dome  of  the  instrument  or  fog  with  visibilities  on  the  order  of  2  km  was 
reported.  R2  increases  to  0.89  if  these  cases  are  removed  from  the  database.  As  indicated, 
the  model  can  not  handle  the  impact  of  fog  on  the  downwelling  infrared  flux.  In  addition, 
surface  temperatures  on  the  order  of  -20  °C  were  reported  in  seventeen  of  the  com¬ 
parisons.  It  is  doubtful  that  the  model  will  predict  the  correct  cloud  base  temperature, 
especially  the  cloud  base  temperature  of  low  clouds,  under  these  conditions.  Usually,  in  a 
winter  environment  when  the  surface  temperature  is  very  low,  there  is  a  strong  inversion. 
Not  only  will  the  calculation  of  the  cloud  base  temperature  be  affected  by  the  presence  of 
an  inversion,  but  the  calculation  of  the  atmospheric  emissivity  will  also  be  affected. 


Table  3.2.3  Grayling  II  average  downwelling  IR  flux  for  different  cloud  conditions. 


Cloud  amount 

Measured  IR 

Calculated  IR 

clear 

203.9 

182.1 

scattered  (0.1  to  0.5  cloud  cover) 

222.1 

225.0 

broken  (0.5  to  0.9  cloud  cover) 

257.5 

269.9 

overcast  (1.0  Cloud  cover) 

289.2 

301.5 

y  =  1.0831X-  17.479 
R2  =  0.8345 
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3.3  Potential  Modification 


As  indicated  above,  there  are  no  provisions  for  handling  inversions  and  optically  thin 
clouds.  With  this  in  mind,  it  is  feasible  to  recast  the  formulation  of  the  cloud  radiance  and 
the  attenuation  of  the  cloud-emitted  radiance  in  terms  of  parameters  that  would  be  readily 
available.  It  should  be  remembered  that  the  proposed  approach  is  only  a  surrogate  for  a 
more  complete  representation  of  cloud  and  atmospheric  processes  that  account  for  the 
downwelling  radiance  at  the  surface.  To  accurately  model  these  more  complicated 
processes,  it  is  necessary  to  have  atmospheric  profiles  of  temperature  and  moisture  plus 
profiles  of  the  cloud  optical  properties,  parameters  that  are  not  readily  available.  The 
proposed  model  is  based  on  the  following  assumptions: 

1 .  Low  and  middle  clouds  are  optically  thick.  This  implies  the  cloud  emissivity  is  one, 
while  the  cloud  reflection  and  transmission  are  zero.  It  also  implies  that  clouds  radiate 
at  the  cloud  base  temperature.  This  assumption  is  not  true  for  ice  clouds  and  some  of 
the  thinner  stratus  type  clouds  found  over  ocean  areas.  Ice  clouds,  because  of  their 
relatively  cold  temperatures,  do  not  contribute  significant  amounts  of  downwelling 
radiance.  For  high  clouds  (ice  clouds),  an  emissivity  of  less  than  one  would  be  used. 

2.  The  lower  atmosphere  contributes  the  greatest  amount  of  clear  sky  radiance  and 
causes  the  greatest  amount  of  attenuation  of  the  cloud-emitted  radiance.  This  implies 
that  an  effective  atmospheric  emissivity  can  be  computed  from  the  ambient  relative 
humidity  (or  vapor  pressure). 

3.  Cloud  overlap,  for  multi-layered  cloud  systems,  is  governed  by  the  principle  of 
random  overlap.  In  reality  the  dynamic  atmospheric  processes  most  likely  govern  the 
overlap.  If  the  dynamic  processes  are  synoptic  scale  in  nature,  they  influence  the 
atmosphere  from  the  surface  to  the  tropopause  and  above.  In  this  case,  the  multi¬ 
layered  clouds  are  most  likely  correlated. 

The  general  equation  for  the  downwelling  radiation  can  be  given  as 

I,i  =  eaoTa  +  (1  - sa)c,aTcl  +  (1  -  sa)c'^Tcm  +  (1  - sa)cheffcrTch  .  (3.13) 

The  effective  atmospheric  emissivity  and  effective  middle  and  high  cloud  cover  are 
calculated  as  indicated  above.  The  cloud  radiating  temperature  is  calculated  as  follows: 

Tc(i,m,h)  =  Ta+T(lat, season, z)Z  (3.14) 

where  T(lat,  season,  z)  is  the  integrated  atmospheric  lapse  rate  and  is  a  function  of 
station  location  (latitude),  season,  and  cloud  base  altitude.  The  lapse  rate  can  be  obtained 
from  the  Geophysics  Laboratory  model  atmospheres  or  from  the  Atmospheric  Circulation 
Statistics  model  atmospheres  of  Oort  and  Rasmusson  (1971).  The  Geophysics  Laboratory 
model  atmospheres  is  stratified  into  three  latitude  zones  and  two  seasons,  while  the 
Atmospheric  Circulation  Statistics  model  atmospheres  is  stratified  by  month  and  five 
degree  latitude  bands.  For  each  latitude  zone  and  each  time  period,  the  integrated  lapse 
rate  could  be  computed  for  the  altitude  range  0-1  km,  0-2  km,  0-3  km,  and  0-10  km.  In 
addition,  the  atmospheric  profile  could  be  adjusted  in  the  lower  levels  of  the  atmosphere 
to  be  consistent  with  the  observed  ambient  temperature  and  humidity.  The  cloud  height, 
if  not  available  from  observations,  can  be  specified  in  a  physically  consistent  manner 
with  the  model  vertical  temperature  profiles. 
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The  impact  of  fog  or  haze  on  the  downwelling  infrared  flux  would  be  parameterized 
using  MODTRAN  (MODerate  resolution  TRANsmission)  and  tied  to  the  visibility  and 
the  cause  of  the  restricted  visibility,  parameters  reported  on  standard  meteorological 
observations.  MODTRAN  is  the  Department  of  Defense’s  standard  atmospheric 
IR/VIS/UV  radiance  and  transmission  band  model  for  lower  altitudes. 


3.4  Upwelling  Longwave  Radiation 

The  upwelling  radiance  can  be  specified  in  terms  of  the  surface  temperature, 
downwelling  radiance,  and  the  surface  emissivity  ( s  ).  If  the  surface  temperature  is  not 

available  it  may  be  necessary  to  assume  the  surface  temperature  is  approximately  equal  to 
the  ambient  temperature.  In  addition,  the  surface  emissivity  is  parameterized  in  tenns  of 
the  land  cover.  The  values  of  the  land  cover  emissivities  can  be  obtained  from  the 
literature.  Assuming  land  cover  (or  land  use)  category  is  available  from  GIS  infonnation. 
The  upwelling  radiance  is  given  as 

(3.15) 


3.5  Model  Input  and  Output 

The  model  inputs  are 

1 .  The  latitude  of  the  location. 

2.  The  Day-Of-the-Year  (DOY)  and  hour  of  the  day  for  the  model  calculations. 

3.  The  land  surface  category.  The  land  surface  categories  currently  modeled  are  grass, 
all  USCS  soil  types,  and  snow.  It  is  assumed  that  the  land  category  is  input  by  the 
user. 

4.  Ambient  (air)  temperature.  The  model  will  not  run  unless  this  parameter  is  provided. 

5.  Relative  Humidity.  If  this  is  missing  the  model  assumes  a  value  of  60%. 

6.  Surface  Temperature.  If  the  surface  temperature  is  missing,  it  is  set  equal  to  the  air 
temperature.  For  daily  average  conditions  this  assumption  is  not  that  bad  since  at 
times  during  the  day  the  surface  temperature  is  greater  than  the  air  temperature  while 
at  other  times  it  is  less  than  the  air  temperature. 

7.  Low,  middle,  and  high  fractional  cloud  amount.  If  all  three  cloud  amounts  are 
missing,  the  low  cloud  amount  is  set  equal  to  the  global  average  cloud  amount  of  0.5. 

8.  Low,  middle,  and  high  cloud  type.  If  the  cloud  type  infonnation  is  missing  and  the 
cloud  amount  is  not  zero,  a  cloud  type  is  assigned  based  on  climatological  averages. 

9.  Low,  middle,  and  high  cloud  height.  If  the  cloud  base  altitude  is  missing  and  the 
cloud  amount  is  not  zero,  a  cloud-base  altitude  is  calculated  based  on  the  cloud  type, 
location  (latitude),  DOY,  and  climatological  infonnation. 

The  model  outputs  are  the  downwelling,  upwelling,  and  net  longwave  radiation  (WZm2)  at 

the  surface  for  the  input  conditions.  If  the  input  conditions  are  representative  of  some 

time  period,  then  the  number  of  Joules/ ni  for  that  period  is  given  as  the  flux  (W/m2) 

times  the  length  of  the  period  in  seconds. 
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Chapter  4 

Soil  Temperature  Model 


4.0  Introduction 

The  soil  temperature  at  the  current  time  step  is  a  function  of  the  past  soil  temperature,  soil 
moisture  state,  and  the  surface  forcing.  The  latter  depends  on  the  time  of  day,  wind, 
relative  humidity,  precipitation,  air  temperature,  and  presence  of  clouds. 


4.1  Soil  Thermal  Model 


The  temperature  gradient  in  a  non-uniform  soil  layer  can  be  described  by  the  one¬ 
dimensional  heat  flow  equation: 
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where  T  is  the  temperature  (K),  t  is  time  (5),  k,h  is  the  thermal  diffusivity  (m  / s ),  cp,w  is  the 
specific  heat  of  water  ( J/kg-K ),  cp  is  the  specific  heat  of  the  soil  ( J/kg-K ),  v  is  the  vertical 
rate  of  water  flow  (m/s),  lfus  is  the  latent  heat  of  fusion  ( J/kg ),  6,  is  the  volumetric  ice 
content  (cm3 /cm3),  p,  is  the  density  of  ice  (kg/m3),  pw  is  the  density  of  water  (kg/m3),  and  z 
is  depth  (m)  measured  positive  downward  from  the  surface.  Further  discussions  on  the 
water  flow  rate  and  the  change  in  ice  content  are  found  in  Chapter  5  and  Chapter  6, 
respectively.  The  temperature  is  subject  to  the  following  boundary  condition: 

F(T)  =  (\-a)Isl+li-l]r  +  H  +  L-P  + 
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(4.2) 


where  a  is  the  surface  albedo,  I s  X  is  the  total  solar  radiation  (W/rn)  impinging  on  the 

1  2  T 

surface,  I.r  is  the  incoming  longwave  radiation  ( W/in  ),  Iir  is  the  outgoing  longwave 

radiation  (W/m  )  both  reflected  and  emitted,  H  is  the  sensible  heat  (W/rn),  L  is  the  latent 
heat  (W/m2),  P  is  the  heat  due  to  precipitation  (W/m2),  and  k  is  the  surface  thermal 
conductivity  (W/m-K).  In  Equation  (4.2),  heat  that  is  transferred  to  the  surface  is 
considered  positive,  as  is  shown  in  Figure  4.1.  Energy  that  can  be  transferred  either 
to/from  the  surface  depending  on  the  gradient  is  shown  as  double -headed  arrows. 


The  first  tenn  in  Equation  (4.2)  represents  the  amount  of  solar,  or  shortwave  radiation, 
absorbed  by  the  surface.  The  second  term  is  the  incoming  longwave  radiation  while  the 
third  term  is  the  outgoing  longwave  radiation,  which  is  composed  of  radiation  emitted 
from  the  surface  and  incoming  radiation  reflected  from  the  surface.  The  sensible  and 
latent  heat  fluxes  together  are  called  the  turbulent  heat  fluxes  and  have  non-zero  values  in 
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the  presence  of  wind.  The  precipitation  heat  represents  the  energy  needed  to  cool  or  heat 
any  snow  or  rain  that  falls  on  the  surface.  The  first  tenn  in  the  second  row  of  Equation 
(4.2)  takes  care  of  heat  conduction  to/from  the  surface  by  the  underlying  ground 
depending  on  the  temperature  gradient.  This  is  followed  by  the  heat  released/absorbed  by 
the  soil  as  the  soil  moisture  melts/freezes.  Finally,  the  last  tenn  represents  heat  that  is 
advected  away  from/towards  the  surface  as  a  result  of  the  vertical  movement  of  moisture. 
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i  A 


\  /><*!,  X 

z  =  o 

H 

L  ‘ 
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.  fus  Pw  dt 

dT 

K  - 

dz  j  . 

Layer  1 
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Layer  3 

Figure  4.1  Surface  Energy  Balance. 


4.1.1  Solar  Radiation 

Both  diffuse  ( Is  Xdif )  and  direct  ( I s  Xd,r )  solar  radiation  comprise  Is  X  such  that 
Is  X=  Is  Xdl'  +Is  Xdi) .  Since  most  surfaces  are  not  horizontal,  the  effect  of  slope  on  the 

amount  of  solar  radiation  that  impacts  a  site  must  be  taken  into  consideration.  The 
position  of  the  sun  is  determined  by  the  solar  zenith  angle,  >pa_  and  the  solar  azimuth 
angle,  i90  as  measured  clockwise  from  North  to  its  horizontal  projection.  Similarly,  the 

slope  is  defined  by  (p,  the  zenith  (elevation)  angle  of  the  surface  as  measured  in  the 
positive  direction  upwards  from  the  horizontal  and  3 ,  the  azimuth  or  aspect  angle 
measured  positive  clockwise  from  North  to  its  horizontal  projection.  Refer  to  Figure  4.2 
for  clarification.  Following  the  work  of  Jordan  (1991)  and  Shapiro  (1978)  the  following 
correction  is  applied  during  daylight  hours,  i.e.,  <p{]  <  90° , 
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I  l  =/  idir 

s  slope  s 


I  Vif 


COS  (p  + 


sin  cp0  sin  cos  Ai9 


+ 


COS  (p0 

(n  -  AS)(Sc+cos<p)  +  A9(\  +  Sccos(p) 
*Q  +  Se) 


+ 


I,  i 


a(l-cos(p ) 


(4.3) 


where  all  angles  are  in  radians  and  A9  =  |i9-i90|  and  Sc  =  1.0  +  0.5 sin (p0  +  2.0sin(2^0)  is 

the  ratio  of  the  average  diffuse  radiance  from  the  solar  and  anti-solar  quadraspheres 
(Jordan  1991).  The  last  term  in  Equation  (4.3)  represents  solar  radiation  reflected  back 
onto  the  surface  due  to  the  slope. 


Figure  4.2  Sloped  Solar  Radiation  Geometry. 


4.1.2  Longwave  Radiation 

All  objects  radiate  energy  over  the  entire  spectrum  proportional  to  their  surface  properties 
and  temperature  according  to  the  Planck’s  function,  l]r  emit  =  scjT 4  where,  s  is  the  broad 

band  surface  emissivity,  T  (K)  is  the  surface  temperature,  and  a  is  the  Stefan-Boltzmann 
constant  (5.669e-08  W/m2 K4).  The  incoming  longwave  radiation  is  either  measured  or 
calculated  according  to  Chapter  3. 

The  outgoing  longwave  radiation  is  composed  of  an  emitted  component,  which  follows 
the  Stefan-Boltzmann  behavior,  and  a  reflected  component,  which  is  proportional  to  the 
incoming  longwave  radiation.  According  to  Kirchoffs  law,  the  emissivity,  transmis¬ 
sivity,  and  reflectance  of  a  surface  must  sum  to  1.0.  The  transmissivity  of  soil  surfaces  is 
approximately  zero,  therefore, 

ll=scjTA  +  {\-s)ll.  (4.4) 
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Unlike  the  shortwave  radiation,  no  corrections  are  made  for  sloped  surfaces. 


4.1.3  Sensible  Heat 


The  sensible  heat  flux  incorporates  the  transfer  of  energy  between  the  surface  and  the 
atmosphere  due  to  molecular  conduction  and  turbulent  mixing.  It  depends  on  the 
temperature  gradient  between  the  two  media.  It  is  given  as  (Jordan  1991) 

#  =  («o+A<>,c»(r,-r)  (4.5) 

with  eo  the  windless  exchange  coefficient  for  sensible  heat  (2.0  W/m  )  (Jordan  1996),  pa 
the  air  density  (kg/m3),  cPyU  the  specific  heat  of  dry  air  at  constant  pressure  (1005.6 
J/kg-K),  Cd  the  dimensionless  drag  coefficient  (discussed  below),  W  the  wind  speed  (m/s) 
at  a  height  of  2.0  m,  and  Ta  the  air  temperature  (K)  at  a  height  of  2.0  m.  Based  on  the 
work  of  Balick  et  al.  (1981),  pa  =  0.00348(f),  /Ta)  where  0.00348  (kg-K/m3-Pa)  is  the 

molecular  weight  of  air  divided  by  the  universal  gas  constant  and  Pa  is  the  air  pressure 
(Pa). 


The  parameterization  of  the  drag  coefficient,  ChD ,  is  dependent  on  the  snow  depth.  If  the 
snow  depth  is  greater  than  measuring  height,  Za  (m)  of  the  air  temperature,  wind  speed 
and  relative  humidity,  then  ChD  =  0.002 +  0.006(Z/ 5000) ,  where  Z  is  the  site  elevation 
(m).  Otherwise  (Koenig  1994), 

C»=ThCf:  (4.6) 

where 


r 


h 


1-0  Rib  <  0.0 

(1.0-16.0J?,,)05 

1.0  Rib  =0.0 

1.0 

1.0-5.0J?/fe  0.0  <  Rib  <  0.2 


(4.7) 


Rib  = 


2gZa{Ta~T) 

(Ta+T)W2 


rg°  —  k 

^ hn 


HZa 


'  ch 


r k  is  the  sensible  heat  exchange  stability  correction  factor,  Cb°  is  the  bulk  transfer 
coefficient  near  the  ground,  is  the  bulk  Richardson  number,  g  is  the  gravitational 
acceleration  (9.81  m/s2),  k  is  von  Karman’s  constant  (0.4),  z()  is  the  ground  roughness 
length  (0.0006  m  for  snow/ice,  0.0001  m  for  pavements,  0.001  m  for  all  soils),  and  rc/t  is 
the  turbulent  Schmidt  number  (0.63). 
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4.1.4  Latent  Heat 


The  latent  heat  flux  term  quantifies  the  energy  lost  or  gained  from  the  surface  due  to 
evaporation  or  condensation,  respectively.  In  order  for  latent  heat  to  be  present,  there 
must  be  a  vapor  gradient  between  the  surface  and  the  atmosphere,  ft  is  formulated  as 
(Kahle  1977,  Balick  et  al.  1981b) 

L  =  paCeDW'l(qg-qa)  (4.8) 

where  pa  is  defined  in  Sec.  4.1.3,  /  is  either  the  latent  heat  of  evaporation,  levap  (2.505e05 
J/kg ),  or  sublimation,  lsub  (2.838e06  J/kg ),  depending  on  the  air  and  surface  temperatures, 
qg  is  the  mixing  ratio  of  the  air  at  the  surface,  and  qa  is  the  mixing  ratio  of  the  air  at  2.0  m, 
and  W'  is  the  corrected  windspeed  (m/s).  If  W  is  below  2.0  m/s  then  W'  =  2.0  m/s  else 
W’  =  W  (Kahle  1977,  Hughes  et  al.  1993).  CeD  is  described  similarly  to  ChD  except  that 


the  turbulent  Schmidt  number,  rcb,  is  replaced  by  the  turbulent  Prandtl  number,  rce  (0.71), 
in  the  expression  for  the  bulk  transfer  coefficient  near  the  ground  and  W1  replaces  W  in 
calculating  the  bulk  Richardson  number  in  Equation  (4.7).  Balick  et  al.  (1981)  give  the 
latent  heat  of  evaporation  as 


Lap  =  2,500, 775.6 -2369.729 


T-T 


-273.15 


(4.9) 


The  mixing  ratio  qa  =  M  qs(T )  with  Mg  the  moisture  factor  ( 0  <  M ,,  <1)  and  qs  the 


saturated  mixing  ratio.  The  value  assigned  to  the  moisture  factor  depends  on  the  degree 
of  saturation  of  the  soil.  If  it  is  raining,  Mg  =  1  otherwise  it  is  equal  to  the  surface  soil 


moisture  content  (refer  to  Chapter  7).  Balick  et  al.  (1981)  quantify  the  mixing  ratios  as 


ds.a  = 


0.622ed 

P-e, 


ed  =RT7  -610.78exp 


17.269(7; -273.15) 
T-  35.86 


(4.10) 


where  is  the  vapor  pressure  (Pa)  and  RHsa  is  the  relative  humidity.  RHS  =  1.0  and 
RHa  =  RH  and 


[  17.269  over  water 
[2 1 .8745  over  ice /  snow 


[35.86  over  water 
[  7.66  over  ice / snow 


4.1.5  Precipitation  Heat 

Jordan  (1991)  quantifies  the  precipitation  heat  flux  as 
P  =  UpcpTp 

fallrate(m  /  hr) 


U  =  —Y 

P  1  P 


3600(sec/  hr) 


(4.11) 
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with  Tp  the  precipitation  temperature  defined  as  the  wet-bulb  temperature,  cp  either  the 
specific  heat  of  water,  cPiW  (4217.7  J/kg-K)  or  ice,  cpj  (-13.3  +  7.807^  J/kg-K)  depending 
on  Tp,  Up  the  mass  precipitation  flux  {kg/m2 -s),  and  y  the  precipitation  density  {kg/m3).  It 
is  assumed  that  the  precipitation  temperature  is  the  same  as  the  air  temperature. 


4.1.6  Final  Top  Boundary  Equation 

Combining  Equations  (4.2),  (4.4),  (4.5)  (4.8)  and  (4.11),  the  top  boundary  condition  is 
reconfigured  as 

F(T)  =  (1  -  a)I,  l  +£/,p  -  eaT'  +  (e,  +  p,cr^W)( T„  -  T) 

-</J  +  ,-/,.c:,7’„  +  k'T7  @  z  =  0m.  (4.12) 

n  r)G 

+lfus-^rte-vcpT  =  0 
Pw  St 

This  is  solved  at  each  time  step  for  the  surface  temperature. 


4.2  Numerical  Solution 

Equation  (4.1)  is  solved  using  a  modified  second-order  Crank-Nicholson  approach. 
Following  the  technique  presented  in  Hombeck  (1975),  Equation  (4.1)  is  rewritten  as 


T  -T 

7+1,1  77  _ 

1 

.  ;s- 

+ 

1 

J5- 

_ 1 

1 

1 

+ 

i _ 

At 

Az 

Az 

y+u+i 


2  T  +T 

7+1,1  T  1  7+l.i- 1 


(M 


+ 


T  -2T  +T 

1 7,1+1 


(Az)2 


(4.13) 


p.w 


7,1 


T  -T 

1 7,1+1  1j,i 

Az 


|  Ifus  Pi 

CP  Pw  At 


where  the  subscripts  j  and  i  represent  time  and  depth,  respectively.  Combining  like  terms 
and  rearranging  so  that  all  terms  involving  7}+ /  are  on  the  left-hand  side  of  the  equation, 
Equation  (4.13)  becomes 


Tj+u- 1 + y  j,Fj+\,i  +  T]+u+ 1  - 
+  Tj  ,i{Pj,i  +  ajj  ~  $  j .,)  +  Tjj+\  (1 _  Pj.,  +  djj)  +  P,j  = 


(4.14) 


where 
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a  -2  -M- 


Cn.w  A Z 


J’’  J’1  r  L. 

I  PU  S, 


//,  “  «/.,  ~  2‘  A/,/  =  2 


*,,=%  ,  nIt=a,;^^e\y 


K,  ’ 


(4.15) 


Equation  (4.14)  is  solved  for  7}+ /  using  a  Newton-Raphson  technique  so  that  the  final 
matrix  equation  becomes 


-yj. i  -1 

-7;, 2  -1 


(4.16) 


where 


-1  -rM- 1  -1  Ar,+i,„ 

5,.  =  -ri+1;M  -  7/  ;  -  r;+y+1  +  ^  2  </</;-  I  . 


7+1,11-1  ^ n-\ 


(4.17) 


7]  must  satisfy  the  boundary  condition  set  forth  in  Equation  (4.2).  Following  the 
technique  of  Kahle  (1977),  Equation  (4.2)  may  be  rewritten  in  the  form 

T4+ClT  +  c2=  0  (4.18) 

with 

c,=—  («„  +  P.c„Ct„w)  +  paC’„W'lMt^t  +U  ,,cp  +  K, 


( 1  -  a)  /,  i  +£/,*  +  (e„  +  /V„C>)  r„  +  upCpr 


a+K\' 


(4.19) 


1  -paCDW'b  Mg  qt(T]^)-^  ■  Tj_u  -qa 

V  V  T=Tj-i,i  J  J 

The  above  assumes  that  kcIT  /  dz  =  k]  (T2  -  7j )  /  Az,  and  that  the  saturated  mixing  ratio 
may  be  represented  as 

g,  =  «,(W+ fjKWfo  -r/-u]  •  H.20> 

Equation  (4.18)  is  solved  for  7  using  the  Newton-Raphson  technique. 

To  ensure  numerical  stability  when  solving  Equation  (4.16),  the  following  criterion  holds 
on  the  time  step  for  each  node: 


At  J_ 

th,i  /  .  \2  ^ 

(Az,)  2 


(4.21) 
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The  thermal  diffusivity  can  change  as  a  function  of  depth  and  time  depending  on  the  soil 
layer  structure.  To  simplify  calculations,  the  vertical  node  increments,  Azt,  are  held  fixed. 
The  time  increment,  At,  is  allowed  to  vary  according  to 

i  At 

ktiui  —  rr 

(Az,) 

Ensuring  that  the  conditions  set  forth  in  Equation  (4.21)  are  met,  a  lower  value  of  0.45  is 
used  to  calculate  the  time  step.  Detennination  of  the  soil  thennal  diffusivity,  kth,  and  other 
physical  properties  needed  to  solve  Equation  (4.16)  are  discussed  in  Section  4.3.  At  each 
calculation  time  step,  At  is  updated  as  necessary  depending  on  the  current  soil  physical 
properties. 

4.3  Material  Properties 

To  determine  the  vertical  temperature  profile  in  the  ground  as  governed  by  Equation 
(4.16),  several  physical  characteristics  of  the  material  are  needed.  They  are  the  surface 
albedo,  a,  and  emissivity,  s,  and  the  nodal  thermal  conductivity,  k,  and  thennal 
diffusivity,  k.  The  latter  two  parameters  are,  among  other  things,  dependent  on  the 
temperature  and  moisture  content  of  the  soil. 

Following  the  recommendations  of  Farouki  (1981),  Johanson’s  method  (1975)  is  used  to 
calculate  the  soil  thermal  conductivity  as  a  function  of  soil  type,  porosity  (n),  dry  density 
(ya),  degree  of  saturation  (S,),  quartz  content  (q),  dry  thermal  conductivity  (tQryX  saturated 
thermal  conductivity  ( Ksat),  thermal  conductivity  of  the  soil  solids  (/y),  unfrozen  water 
volume  {ww),  and  the  Kersten  number  (Ke).  The  governing  equation  is 

K  ( Ksat  —  K(fry)Ke  Kdty  (4.23) 

where 

provided  by  user 
13»'  +  64-7'  ±M% 

2700-947yrf  (4.24) 

0.05  Tj  .  >  273. 15 A  peat 

0.55  Tj,  <273.1 5K  peat 

0.57"k*(1~")  T..  >273.1 5K 

A  J  d 

2.2"/c(1“")0.269ww  T. .  <  273. 15A 

'  ]’1  (4.25) 

0.55  Tj  t  >273.15^  peat 

1.80  Tjj  <273.15A  peat 

0.7 log ^  +  1.0  Tjj  >  273. 15A  coarse 

Ke  =  \  log  Sr  +  1 .0  Tjj  >  273 . 1 5 K  fine  (4.26) 

5,  Tjj  <  273. 15 AT 
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K 


[7.7*3. 0(1_,)  q  <  0.20  coarse  .  . 

{  .  (4.27) 

[7.7?2.0(1_9)  a//  others 

and  rM,  =0.001 877 y  +0.0787,  /r;.  =488.197^ +0.4685  .  Farouki  (1981)  found  that 

Johansen’s  method  for  predicting  soil  thennal  conductivity  matched  the  existing  data 
well  except  for  the  case  of  dry  (Sr  <  0.2),  unfrozen,  coarse  soils.  Since  the  recommended 
methods  for  this  situation  rely  on  unavailable  curve  fitting  parameters,  Johansen  (1975)  is 
used  for  all  conditions.  The  physical  properties  used  in  Equations  (4.23)-(4.27)  for  the  15 
USCS  soil  types  are  given  in  Table  4.3.1. 


The  soil  thermal  diffusivity,  k,h,  is  equal  to 

**=§  (4-28) 

where  the  volumetric  heat  capacity,  C,  is  defined  as  (Farouki  1981) 

C  =  (l- n)rdCs  +{ww) pwCw  +{n - [ww  +  wi])paCa  +  (wi)pi  (-13.3  +  7.87V )  (4.29) 

and  Cw  =  4217.7  J/kg-K  is  the  specific  heat  of  water,  Ca  =  1250.0  J/kg-K  is  the  specific 
heat  of  air,  and  Cs  is  the  specific  heat  of  solids  given  in  Table  4.3.1. 


Table  4.3.1  De 


fault  Soil  Properties 


Soil  Type 

Bulk 

Dry  Density 

yd  (g/cm3) 

Porosity 

n 

Emissivity 1 
£ 

Albedo 1 

a 

Quartz 
Content 2 

q 

Organic 

Fraction 

0of 

Specific 

Heat3 

Cs  ( J/kg-K) 

Fine/Coarse* 

GW 

1.95s 

0.296s 

0.92 

0.40 

0.65 

0.00 

820.0 

1 

GP 

2.16s 

0.203s 

0.92 

0.40 

0.65 

0.00 

820.0 

1 

GM 

1.91s 

0.324s 

0.95 

0.40 

0.65 

0.00 

820.0 

1 

GC 

1.87s 

0.344 

0.92 

0.40 

0.65 

0.00 

820.0 

1 

sw 

1.876s 

0.320s 

0.92 

0.40 

0.80 

0.00 

830.0 

1 

SP 

1.594s 

0.415s 

0.92 

0.35 

0.80 

0.00 

816.47 

1 

SM 

1.474s 

0.526s 

0.92 

0.35 

0.80 

0.05 

850.6 

1 

SC 

1.88s 

0.400s 

0.92 

0.35 

0.80 

0.05 

830.0 

1 

ML 

1.457s 

0.464s 

0.94 

0.40 

0.35 

0.10 

845.7' 

2 

CL 

1.589s 

0.422s 

0.97 

0.23 

0.05 

0.10 

854.1 7 

2 

OL 

1.165s 

0.533s 

0.955 

0.265 

0.20 

0.25 

837.4' 

2 

CH 

1.517s 

0.457s 

0.98 

0.30 

0.05 

0.10 

845.7' 

2 

MH 

1.0601 

0.5471 

0.94 

0.30 

0.35 

0.10 

830.0 

2 

OH 

0.841 1 

0.8924 

0.955 

0.265 

0.20 

0.25 

866.7' 

2 

PT 

0.25 

0.70 

0.92 

0.40 

0.05 

0.50 

830.0 

2 

SMSC 

(MC) 

1.601 

0.3961 

0.92 

0.40 

0.80 

0.05 

830.0 

1 

CLML 

(CM) 

1.617s 

0.397s 

0.96 

0.30 

0.20 

0.10 

830.0 

2 

EV 

1.876s 

0.320s 

0.92 

0.40 

0.80 

0.00 

830.0 

1 

CO 

2.185 

0.020 

0.906 

0.406 

0.00 

850.0 

AS 

2.500 

0.020 

0.946 

0.1256 

0.00 

880.0 

RO 

2.700 

0.020 

0.89 

0.40 

0.00 

800.0 

SN 

0.920 

0.020 

0.90 

0.70 

0.00 

1  Sullivan  et  al.  (1997),  p.  28,  38,  and  51 

2  Tarnawski  et  al.  (1997),  p.  96 
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3  Farouki  (1981),  p.  136 

4  Yuet  al.  (1993) 

5  Guymon,  et  al.  (1993),  p.  21,  25,  47-59 

6  Balick,  L.K.  et  al.  (1981),  Table  B1 


4.4  Model  Validation 

We  conducted  a  comparison  of  the  surface  temperatures  generated  using  the  model  and 
measured  values  for  Grayling,  MI,  during  the  fall  of  1992  and  for  Yuma,  AZ,  during  the 
spring  of  1993  and  Hanover,  NH,  during  the  winter  of  2002.  The  data  collected  consist  of 
soil  temperature  profiles  as  well  as  meteorological  conditions.  Almost  no  soil  properties 
data  were  collected,  thus  we  had  to  estimate  these  based  on  values  found  in  the  literature 
for  the  general  soil  types  at  the  sites.  The  Yuma  and  Grayling  locations  are  sandy  while 
the  CRREL  site  is  pavement.  The  data  collected  for  the  first  two  sites  are  part  of  the 
SWOE  series. 


Yuma,  AZ 

|  probe  1  probe  2  ^"FASST  | 


Figure  4.3  Comparison  of  model  results  and  measurements  for  surface  temperature  (K) 
for  Yuma,  AZ. 

Figure  4.3  demonstrates  the  capability  of  FASST.  Except  for  soil  type,  no  other 
information  was  collected.  The  two  temperature  probes  were  placed  side  by  side  in  the 
vicinity  of  the  meteorological  station. 

The  initial  soil  temperature  profile  is  either  input  by  the  user  or  generated  by  the  model, 
which  assumes  a  uniform  profile  equal  to  the  air  temperature.  In  either  case,  the  model 
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assumes  that  the  bottom  heat  flux  is  constant  throughout  the  exercise.  One  example  of  the 
difference  between  having  measured  initial  soil  temperature  profile  data  and  letting  the 
program  calculate  one  is  seen  in  Figure  4.4  below  for  the  Yuma,  AZ,  data.  The  “yes”  case 
is  the  same  as  the  model  results  shown  in  Figure  4.3.  Based  on  the  user-specified  number 
of  layers,  position  and  number  of  initial  soil  temperature  and  moisture  measurements,  and 
the  total  thickness,  the  model  divided  the  profile  into  15  nodes.  For  the  two  “no”  cases, 
everything  was  retained  from  the  “yes”  scenario  except  the  measured  initial  soil 
temperatures.  In  one  “no”  run,  the  number  of  nodes  was  forced  to  be  the  same  as  in  the 
“yes”  case;  in  the  other  “no”  run,  the  model  determined  the  required  number.  As  can  be 
seen,  essentially  no  difference  is  discerned  between  the  three  runs  except  at  the  start  of 
the  simulation. 


Yuma,  AZ  1993 

|™yes,  15  ^Tio,  15  ^Tio,  12 
330 


325 


285 


280  4 - T - T - 1 - T - 1 - T - 1 - T - 1 - 1 

73  74  75  76  77  78  79  80  81  82  83 

Day  of  the  Year 

Figure  4.4  Effect  on  the  results  of  having  initial  soil  temperature  profile  data  (yes)  versus 
letting  the  model  compute  one  (no).  The  number  in  the  legend  indicates  the  number  of 
nodes. 

We  used  the  Yuma,  AZ,  data  to  perform  a  sensitivity  study  on  the  model.  The  variables 
tested  were  number  of  layers,  albedo,  and  emissivity.  We  let  the  model  calculate  the 
thermal  conductivity  and  diffusivity.  We  found  no  dependence  in  the  results  on  the 
number  of  layers.  The  results  showed  little  sensitivity  to  the  emissivity.  The  albedo  does 
affect  the  surface  temperature  since  the  larger  the  albedo,  the  less  shortwave  radiation  is 
absorbed  by  the  surface.  The  results  are  seen  in  Figure  4.5. 
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Yuma,  AZ  1993  Albedo  Test 


1  —  0.2  0.3  0.4  —  0.5 


Figure  4.5  Sensitivity  of  the  surface  temperature  to  the  albedo. 


Grayling,  Ml  1 


- SW  - SM  - SC  - SP 


Day  of  the  Year 

Figure  4.6  The  effect  of  soil  type  on  surface  temperature  for  Grayling,  MI. 
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Since  the  category  of  “sand”  is  broad,  we  investigated  the  effects  of  using  different 
combinations  of  sand  layers  with  the  Grayling,  MI,  data.  We  tested  all  four  types:  SM 
(silty  sand),  SW  (well-graded  sand),  SP  (poorly  graded  sand),  and  SC  (clayey  sand).  The 
results  are  presented  in  Figure  4.6.  The  differences  follow  the  thermal  conductivity 
gradient,  with  the  SW  soil  having  the  lowest  thermal  conductivity  and  the  SM  soil  the 
highest. 

FASST  will  also  be  used  to  predict  the  surface  condition  of  pavement  for  mobility 
purposes.  As  can  be  seen  in  Figure  4.7,  it  does  this  very  well  for  the  data  collected  at  the 
Army  Corps  of  Engineer’s  CRREL,  Hanover,  NH,  site,  especially  when  compared  to  the 
spot  radiometer,  which  is  more  accurate  than  the  thermistor. 


CRREL  Pavement  Study,  2002 

[^"measured  ^"FASST  ^"spot  radiometer  | 


300.0  n 


54  55  56  57  58 

Day  of  the  Year 

Figure  4.7  Comparison  of  FASST  and  measured  pavement  surface  temperatures. 
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Chapter  5 

Soil  Moisture  and  Strength  Models 


5.0  Introduction 

Calculating  changes  in  the  volumetric  moisture  content  profile  as  a  function  of  time  is  the 
critical  step  in  determining  the  state  of  the  ground.  This  is  in  part  due  to  the  highly 
nonlinear  nature  of  the  parameters  associated  with  soil  moisture. 


5.1  Soil  Moisture 


5.1.1  Governing  Equations 


The  flow  of  water  (v)  through  a  porous  media  is  governed  by  Darcy’s  Law,  which  states 
that 


Ydk 

v  =  K  — 
dz 


(5.1) 


where  K  (m/s)  is  the  hydraulic  conductivity,  z  (m)  is  the  depth,  positive  downward  from 
the  surface,  and  h  (in)  the  total  head  equals  the  elevation  head,  or  depth,  (z)  minus  the 

3 

pressure  head  (i//),  i.e.,  h  =  zcos^-'F  =  zcos<p-Pf  / pwg  ,  Pf(Pa)  is  pressure,  pw  (kg/m  ) 

is  the  density  of  water,  cp  is  the  surface  slope,  and  g  (rn/s)  is  gravity.  If  the  soil  is 
unsaturated,  y/<  0  and  h  >  zcoscp.  For  saturated  soils  i//>  0,  requiring  that  h  <  zt.cos(p. 

Also  governing  the  flow  of  moisture  through  a  soil  is  the  conservation  of  mass,  which 
states  that  the  time  rate  of  change  of  the  moisture  content  in  a  given  volume  equals  the 
net  gain/loss  of  fluid  in  the  volume,  i.e., 

5^  =  _5v_  p.  d6l 
dz 


-  +  sources  -  losses 


(5.2) 


dt  dz  pw  dt 

where  0W  (cm3 /cm3)  is  the  volumetric  moisture  content,  0,  (cm3 /cm3)  is  the  volumetric  ice 

3 

content,  p,  (kg/m)  is  the  density  of  ice,  and  t  (sec)  is  time.  Equation  (5.2)  assumes  that 
changes  with  respect  to  time  in  the  soil  porosity  and  water  density  are  negligible 
compared  to  changes  in  the  soil  moisture  and  total  head.  Further  discussion  on  the  change 
in  ice  content  is  found  in  Chapter  6.  The  source  and  loss  terms  in  this  equation  account 
for  occurrences  such  as  runoff  and  plant  root  uptake.  The  latter  is  for  future  development. 

Equation  (5.2)  is  subject  to  the  following  flow  boundary  conditions  at  the  surface  and  at 
the  bottom  of  the  soil  column: 


q<oP  --E  +  Cr+  P  +  (hpond  +  himelt  +  hsmelt ) / At  @  z  -  0 

qho,  =  K  sin  (slope)  @  z  =  zhot 


(5.3) 
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where  E  (m/s)  is  the  evaporation  rate,  C,  (m/s)  is  the  condensation  rate,  P  (m/s)  is  the  rate 
of  precipitation,  and  hpond  (m)  is  the  head  due  to  water  collecting  on  the  surface,  and  hiymeit 
(m)  and  hsmeit  (m)  are  the  heads  due  to  melting  ice  and  snow,  respectively,  and  At  (sec)  is 
the  time  step.  If  the  ground  is  sloped,  no  water  accumulates  and  any  water  that  falls  on 
the  surface,  but  which  does  not  infiltrate,  becomes  runoff. 


The  evaporation  and  condensation  rates  depend  on  the  moisture  difference  between  the 
air  and  the  surface.  They  are  quantified  as  the  latent  heat  flux.  Following  the  procedure 
outlined  in  Chapter  4,  Section  4.1.6,  the  modified  surface  energy  balance  is  rearranged 
such  that 


(l-a)/,  4.  +sl* - ecTA  +  pacpaCDW iTa - T ) 


latent  heat  = 


(5.4) 


Reference  should  be  made  to  Chapter  6  for  explanation  of  the  individual  tenns.  The 
moisture  flux  to/from  the  surface  is  then 

latent  heat 


moisture  flux  = 


Pj 


(5.5) 


and  /  (J/kg)  is  the  latent  heat  of  evaporation  if  the  air  temperature  is  above  freezing  and 
the  latent  heat  of  sublimation  otherwise.  If  the  latent  heat  flux  is  positive,  evaporation  is 
occurring,  otherwise  condensation  is  happening. 


5.1.2  Numerical  Solution 
Equation  (5.2)  is 


At 

where 

=  Ki- m 

Vi+ 1  =  Kj+1/2  I  I 

L  Zi+\  ~  Zi  J 

hi  =  zi  cos (p -  IT  ,  Azi  =  (zM  -zm)/2,  n,  is  the  porosity  at  i  and  the  subscripts  j  and  i 

represent  time  and  depth,  respectively.  The  change  in  soil  moisture  content  due  to 
changes  in  the  ice  content,  i.e.,  freezing/thawing,  is  incorporated  into  the  source  and  sink 
tenns.  In  Equation  (5.2)  it  is  the  second  term  on  the  right-hand  side.  Equation  (5.6)  is 
solved  for  i//,  using  a  Newton-Raphson  technique  so  that  the  final  matrix  equation 
becomes 


Dived  numerically  using  an  explicit  scheme  such  that 
+  sources(i)  -  losses (i) 


Vj+ l,i+l  Vj+l,i 


A z. 


(5.6) 


i 


zi~zi- 1 


hM  ~hi 


(5.7) 
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(5.8) 


1 - 

+  ^ 

l _ 

rA'F 

7+1,1 

Bl  " 

a2  y  2  u2 

A'F 

^  X  7+1,2 

b2 

a  ,  y  ,  a+ , 

n- 1  /  n- 1  n- 1 

A'F 

^  T  y+l,H-l 

B, 

a~  y 

|_  n  in  J 

A'F 

^  T  7+l.H 

A  _ 

where 


a,  =  ■ 


8B 


W(+ 

SB 


vp  _vp 
cos  cp - — - 


Tf  -f. 


a.  =  ■ 


SHV  ,  Az 


8B 
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1/2 
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’  B>  =  ^  +^r{vj+i,M  -  v/+i,, )  - sourcesii)  +  losses{i) .  (5.11) 


(5.9) 


/  =  1 


2<i<  « -1  (5.10) 


/  =  n 


To  ensure  numerical  stability  when  solving  Equation  (5.8),  At  is  chosen  such  that 
(Ksat  )  At  <  Az; ,  where  (Ksat)j  is  the  saturated  hydraulic  conductivity.  Default  values  for 
Ksat  for  the  different  USCS  soil  types  are  found  in  Table  5.1.1  located  in  the  next  section. 


5.1.3  Hydraulic  Parameters 


The  relationship  between  volumetric  moisture  content  and  pressure  head  is  highly 
nonlinear.  Following  the  work  of  van  Genuchten  (1980)  it  is 

0w  =  0r  + - ^"Iax  ~  6‘ -  (5.12) 

(i+MT 


where  6r  is  the  residual  volumetric  water  content,  6max  is  the  maximum  volumetric  water 
content,  aVG  (cmT1)  is  a  constant  related  to  the  reciprocal  of  the  bubbling  pressure  head,  n 
is  a  constant  dependent  on  the  distribution  of  pores,  and  mvc  =  1  -  1 /«  ,<,'•  Default  values 
for  6r ,  Omax,  ccvg,  and  nVG  are  found  in  Table  5.1.1.  As  can  be  seen  in  Figure  5.1,  the 
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relationship  between  the  volumetric  moisture  content  and  the  pressure  head  is  highly 
nonlinear. 


Table  5.1.1 

Default  Soil  Moisture  Properties. 

Soil  Type 

Saturated  Hydraulic 
Conductivity 

Ksat 

(cm/s) 

Residual 

Water 

Content 

0r  ; 

(m3/m3) 

Maximum 
Water  Content 

0  i !  i ;  l  \ 

(m3/m3) 

Van 

Genucht. exponent 

nvG 

Bubbling 
Pressure  Head 

CCvG 

(cm) 

GW 

0.000262 

0.01 

0.2962 

1.53 

22.61 

GP 

0.00010562 

0.01 

0.4032 

2.23 

22.61 

GM 

0.0000672 

0.01 

0.3242 

1.23 

32. 71 

GC 

0.000013891 

0.01 

0.344 

1.53 

23.21 

sw 

0.000023612 

0.01 

0.3202 

1.253 

38.547J 

SP 

0.000007412 

0.01 

0.4152 

2.53 

38.547J 

SM 

7.987e-062 

0.01 

0.5262 

1.43 

68.5473 

SC 

0.000000292 

0.01 

0.4002 

1.53 

58.5473 

ML 

0.000000572 

0.01 

0.4642 

1.53 

33.91 

CL 

7.7e-082 

0.01 

0.4222 

1.343 

53. 51 

OL 

0.000097222 

0.01 

0.5332 

1.343 

28. 61 

CH 

4.8e-082 

0.01 

0.4572 

1.53 

32.91 

MH 

1 .5e-08 

0.01 

0.5472 

1.343 

39. 01 

OH 

0.0008831 

0.01 

0.8924 

1.343 

29.31 

PT 

0.000014 

0.15 

0.70 

1.34 

38.6 

SMSC 

(MC) 

7.7e-072 

0.01 

0.3962 

1.53 

23.5473 

CLML 

(CM) 

1.26e-072 

0.01 

0.3972 

1.343 

32.91 

EV 

0.000023612 

0.01 

0.3202 

1.253 

38.547J 

CO,  AS, 
RO,  SN 

0.0 

0.001 

0.02 

1  Sullivan  et  al.  (1997),  p.  28 

2  Guymon  et  al.  (1993),  p.  21,  25,  47-59 

3  Jordan  (2000) 

4  Yuet  al.  (1993) 
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K  =  Ksatw' 


0.5 


l-(l -w1^)" 


~l2 


where 


w  = 


6  -6 

w  r 

e . -o, 


(5.13) 


(5.14) 


A  plot  of  Equation  (5.13)  as  a  function  of  0  is  presented  in  Figure  7.2  along  with  the 
hydraulic  gradient,  dQw/dy/ used  in  Equation  (5.11). 


In  a  layered  system,  the  average  hydraulic  conductivities  A)-//?  and  Ki+i/2  found  in 
Equations  (5.7),  (5.9)  and  (5.10)  equal 


K,- 1/2 
Ki+ 1/2 


KiAzi_l+Ki_lAzi 

K:KmAz,Azm  ' 

KM+i+KMAzt 


(5.15) 


Their  derivatives  with  respect  to  the  pressure  head  are  of  the  general  form 


dKj_  1/2  _  dAM/2  GKi  ,  gw,,,  d0Wi  i 

aAf._i/2  _  a^,._1/2  f/c  aw,, 

dTV  5/f,  86w  &¥. 

dKi+ 1/2  _  d^f+l/2  d^f+1  dwf+1  ^ 
^!+1  dKM  dwM  80Wm  dVM 

dKi+ 1/2  _  dif/+1/2  gif,  gw,.  80n, 

stv  a/r,.  a6»H,  shf 


(5.16) 


The  problem  is  now  fully  developed  and  comparisons  may  be  made  between  the  model 
predictions  and  actual  field  measurements. 


5.1.4  Model  Validation 

As  for  the  soil  temperature  model  verification,  we  used  data  from  the  SWOE  program  to 
validate  the  model.  The  specific  locations  are  Grayling,  MI,  and  Yuma,  AZ.  In  both 
cases,  soil  moisture  by  percent  weight  measurements  were  taken  essentially  once  a  day. 
Researchers  at  CRREL  (Cold  Regions  Research  and  Engineering  Laboratory)  and  WES 
(Waterways  Experiment  Station)  made  independent  measurements.  The  results  for 
Grayling,  MI,  are  shown  in  Figure  5.3.  The  percent  weight  measurements  were  converted 
to  volumetric  soil  moisture  values  to  be  consistent  with  the  model  output.  In  doing  so,  we 
assumed  a  soil  density  of  1.49  g/cm3.  Also  for  the  Grayling  comparisons,  an  initial 
surface  soil  moisture  of  0.05  m  /m  was  assumed  based  on  the  measured  data.  The  model 
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calculations  are  well  within  the  spread  of  the  measurements  and  respond  to  precipitation 
events  as  seen  in  Figure  5.3.  Soil  parameters  were  chosen  to  obtain  the  best  fit. 

Grayling,  Ml  1 


♦  WES  1  WES  2  A  CRREL - FASST  - rain  snow 


Figure  5.3  Soil  moisture  comparisons  for  Grayling,  MI,  1992. 

We  used  the  CRREL  daily  measurement  for  Yuma,  AZ,  taken  on  day  74,  the  first  day  of 
the  experiment,  to  initialize  the  surface  soil  moisture  content.  The  moisture  content  at 
depth  was  initialized  using  the  procedure  discussed  in  Chapter  11,  Section  1 1.2.  The  soil 
is  user-specified  SM,  or  silty  sand.  Overall,  the  model  does  fairly  well  predicting  the  soil 
moisture,  as  can  be  seen  in  Figure  5.4.  It  is  certainly  within  the  measurement  error  of  the 
probes  used  by  WES  and  CRREL.  The  daily  fluctuations  are  due  to  changes  in 
evaporation/condensation. 

Since  no  specific  soil  type  except  for  sand  was  recorded  for  Grayling,  we  investigated  the 
effect  of  soil  type  on  the  moisture  content.  As  with  the  surface  temperature  data,  there 
were  only  small  differences  in  the  model  results  for  the  different  soil  types  as  seen  in 
Figure  5.5.  Unlike  the  saturated  hydraulic  conductivity  values  used  to  generate  the  results 
shown  in  Figures  5.3  and  5.4,  the  default  values  used  to  generate  the  results  in  Figure  5.5 
are  much  smaller  and,  therefore,  no  daily  oscillations  are  observed.  In  all  the  figures,  it  is 
evident  that  soil  moisture  measurements  are  inexact,  with  as  much  as  a  10%  difference 
exhibited  on  a  given  day. 
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Figure  5.4  Soil  moisture  comparisons  for  Yuma,  AZ. 


105  110  115  120 


Grayling,  Ml  1 


♦  WES  Speedy  ■  WESTroxel  a  CRREL  SM  SW  SC - SP  - rain  - snow] 


Day  of  the  Year 

Figure  5.5  Investigation  of  soil  type  on  soil  moisture  for  Grayling,  MI,  1992. 
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5.2  Soil  Strength 

The  soil  strength  module  was  developed  for  the  Soil  Moisture  Strength  Prediction  Model 
Version  II  (SMSP  II).  The  original  documentation  is  found  in  Sullivan  et  al.  (1997). 
Section  7.2.1  is  a  paraphrase  of  pages  32-33  and  75-76  of  the  aforementioned  report. 


5.2.1  Theory 

Historically,  cone  index  (Cl)  as  measured  with  a  standard  WES  (USA  Corps  of  Engineers 
Waterways  Experimental  Station)  cone  penetrometer,  has  been  used  to  predict  the 
strength  of  a  soil  in  establishing  vehicle-terrain  interaction  relationships.  The  cone  index 
is  a  measure  of  the  resistance  to  penetration  developed  by  the  cone  as  it  is  pushed  into  the 
soil.  It  equals  the  vertical  force  applied  to  the  sleeve  divided  by  its  surface  area.  The  WES 
cone  has  a  30°  apex  angle  with  a  0.5 -in.  base  area  that  is  connected  to  a  circular  shaft  of 
3/8 -in.  diameter. 

The  rating  cone  index  (RCI),  on  the  other  hand,  has  been  used  to  represent  the  proportion 
of  the  original  soil  strength  retained  after  vehicles  have  passed  over  the  area.  The  rating 
cone  index  is  the  cone  index  multiplied  by  the  remolding  index  (RI)  of  the  soil  where  the 
RI  equals  the  ratio  of  the  after-to-before  Cl  values  of  a  soil  subjected  to  a  remolding  test. 
The  first  step  of  the  remolding  test  is  to  obtain  an  undisturbed  soil  sample  using  a 
trafficability  sampler.  This  sample  is  then  placed  in  a  cylinder  where  the  cone 
penetrometer  tests  are  performed.  Measurements  are  taken  at  1-in.  intervals.  After  placing 
the  soil  in  the  cylinder,  the  sample  is  subjected  to  100  blows  with  a  2.5 -lb  hammer,  12-in. 
fall  per  blow.  A  second  cone  penetrometer  test  is  performed  and  the  resulting  RI  is 
defined  as  the  ratio  of  the  average  Cl  value  after  the  100  blows  to  that  of  the  Cl  value 
before  the  blows. 

Most  soils  have  a  minimum  strength  that  will  permit  a  vehicle  to  complete  a  specified 
number  of  passes  over  it.  This  value  is  referred  to  as  the  vehicle  cone  index  (VCI). 
Insufficient  soil  strength  may  cause  wheeled  and  tracked  vehicles  to  become  immobilized 
if  the  in-situ  Cl  or  RCI  is  less  than  the  VCI  required  for  a  specified  number  of  vehicle 
passes. 

Soil  strength  depends  not  only  on  the  soil  type,  but  also  on  the  soil  moisture.  To  help 
quantify  the  relation  between  soil  strength  and  moisture,  over  1000  field  and  laboratory 
tests  were  perfonned  on  samples  collected  from  the  United  States  and  elsewhere.  Except 
for  gravels  and  peat,  the  resulting  curve  fits  for  the  different  USCS  soils  are 

Cl  =  exp[Cj  -  c2  ln(MC)]  i 

RCI  =  expffj  -c2  ln(MC)] 

where  a  and  c?  are  given  in  Table  7.2.1  below  and  MC  (%)  is  the  soil  moisture  content 
by  weight,  MC  =  0{ pw  /  p. )  •  1 00% ,  and  ps  is  the  soil  density.  The  cone  index  and  rating 

cone  index  for  gravels  is  300,  the  maximum  soil  strength  allowed.  Peat  is  assumed  to 
have  no  strength. 
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The  above  analysis  of  soil  strength  assumes  that  the  soil  is  not  frozen.  Cold  regions  soil 
conditions  present  unique  difficulties  in  characterizing  the  trafficability  of  soil.  A  cone 
index  for  a  soil  with  a  surface  frost  of  1-2  cm  is  irrelevant  because  the  high  Cl  obtained  is 
indicative  of  the  difficulty  in  penetrating  the  frozen  soil  with  a  handheld  instrument.  It  is 
not  necessarily  a  measure  of  the  bearing  capacity  of  the  soil  for  vehicles  or  repeated  foot 
traffic.  Similarly,  the  trafficablility  of  a  soil  with  a  thawed  layer  is  dependent  as  much  on 
thaw  depth  as  it  is  on  the  Cl  of  the  thawed  soil.  A  thaw  layer  so  shallow  that  personnel  or 
vehicles  sink  through  it  to  the  underlying  supportive  frozen  ground  need  not  be 
detrimental  to  mobility,  despite  the  fact  that  the  surface  soil  has  a  low  CL  Another 
complication  is  that  because  an  average  Cl  for  the  top  15  cm  of  soil  is  standard  in 
mobility  applications,  the  bearing  capacity  of  a  thaw  layer  less  than  15  cm  deep  will  be 
misrepresented  by  the  standard  CL 


Table  5.2.1  Defau 

It  Soil  Strength  Coefficients. 

Soil  Type 

Cl  Coeff. 

Cl 

Cl  Coeff. 

C2 

RCI  Coeff. 

Cl 

RCI  Coeff. 

Cl 

GW 

GP 

GM 

GC 

SW 

3.987 

0.81500 

3.987 

0.8150 

SP 

3.987 

0.81500 

3.987 

0.8150 

SM 

8.749 

-1.1949 

12.542 

-2.955 

SC 

9.056 

-1.3566 

12.542 

-2.955 

ML 

10.225 

-1.565 

11.936 

-2.407 

CL 

10.998 

-1.848 

15.506 

-3.530 

OL 

10.977 

-1.754 

17.399 

-3.584 

CH 

13.816 

-5.583 

13.686 

-2.705 

MH 

12.321 

-2.044 

23.641 

-5.191 

OH 

13.046 

-2.172 

12.189 

-1.942 

PT 

SMSC 

(MC) 

9.056 

-1.3566 

12.542 

-2.955 

CLML 

(CM) 

9.454 

-1.3850 

14.236 

-3.137 

EV 

3.987 

0.81500 

3.987 

0.8150 

CO,  AS,  RO,  SN 

Sullivan  et  al.  (1997),  p.  33  and  34 


Despite  these  difficulties,  infonnation  highly  relevant  to  trafficability  can  be  obtained 
from  FASST  based  on  the  predicted  frost  and  thaw  depths  and  the  volume  of  meltwater 
produced  by  the  snow.  The  presence  of  5  cm  of  frost  at  the  surface  will  usually  allow 
unlimited  cross-country  operation  if,  when  unfrozen,  the  soil  is  not  extremely  difficult  to 
traverse  (Richmond  1991). Therefore,  when  FASST  predicts  a  frost  depth  of  at  least  5  cm 
in  a  soil  that  is  predominantly  silt  or  sand,  the  associated  pseudo-CI  is  300,  the  maximum 
strength  allowed.  If  the  terrain  is  normally  untrafficable,  then  a  frost  depth  of  50  cm  is 
needed  in  order  that  most  vehicles  do  not  break  through  the  frost  layer  (Richmond  et  al. 
1995).  For  clays  and  peat  with  predicted  frost  depths  of  at  least  50  cm,  the  pseudo-CI  is 
300. 
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A  finer  resolution  of  bearing  capacity  when  frozen  soil  overlies  incompetent  soil  has  been 
formulated  by  Shoop  (1995).  The  minimum  frost  depth  to  prevent  breakthrough  is 
expressed  in  terms  of  the  vehicle  class  (for  example,  1  person  =  0.1,  pickup  truck  =  4, 
5-ton  truck  =  12,  and  a  D7H  tractor  =  25  [Shoop  1995  and  http://155.217.58.58/cgi- 
bin/atdl.dll/fm/3034.343/appb.htm1).  which  is  the  military  load  classification  for  bridge 
and  highway  limits.  It  can  be  approximated  using  the  gross  weight  of  the  vehicle  (or 
group  of  people)  in  tons.  For  wet  soil  beneath  the  frozen  layer,  the  relation  is 


zwet  =  OAO^vehicle  class 
and  for  dry  soil  the  relationship  is 


(5.18) 


(5.19) 


where  z  (m)  is  the  frost  depth. 

Another  contribution  of  FASST  to  trafficability  is  a  quantification  of  the  amount  of  water 
available  to  infiltrate  thawed  soil  from  a  melting  snow  cover.  Once  the  wetness  of  the 
thawed  layer  is  known,  FASST  can  assess  the  bearing  capacity  of  the  thawed  layer  as  if  it 
were  an  infinitely  thick  layer  (ignoring  the  underlying  frozen  soil).  If  the  bearing  capacity 
of  the  soil  that  forms  the  thawed  layer  is  high,  then  the  presence  of  the  thawed  layer  is 
largely  irrelevant.  If  the  bearing  capacity  of  the  thawed  layer  is  low  or  moderate,  then  its 
thickness  must  be  considered  in  determining  its  impact  on  trafficability. 

Finally,  FASST  outputs  a  “slippery  factor”  that  tells  the  mobility  model  whether  the 
surface  is  dry,  wet,  snow  covered,  or  ice  covered. 
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Chapter  6 

Ground  State  Model 


6.0  Introduction 

The  goal  of  the  Ground  State  Module  is  to  provide  a  very  simple  formulation  and  very 
rapid  calculation  of  freezing  and  thawing  that  still  accounts  for  fundamental  mechanisms 
and  conditions  in  individual  cases.  The  overall  strategy  is  to  use  simple  flux  balance 
relations  to  evaluate  the  interface  temperatures  in  terms  of  energy  input  and  medium 
properties;  to  then  use  those  temperatures  to  express  the  fluxes;  and  then  use  fluxes  at 
depth  to  express  rate  of  freeze  or  thaw.  This  rate,  multiplied  by  the  time  step,  then  gives 
the  depth  of  freeze  or  thaw.  This  depth  of  freeze  or  thaw  can  be  continuous  or  discon¬ 
tinuous,  depending  on  the  temperature  and  moisture  profiles  of  the  soil. 


6.1  Physical  Setup 

The  phase  change  process  is  assumed  to  be  isothermal.  The  main  principle  behind  the 
freeze-thaw  model  is  to  compare  the  energy  needed  to  freeze  or  thaw  the  soil  to  that 
available  at  a  given  depth  as  was  done  by  Guyman  et  al.  (1993).  Partial  freezing  or 
thawing  can  occur.  If  the  ice  content  at  a  given  node  is  greater  than  zero,  the  node  is 
assumed  frozen  for  the  purposes  of  calculating  the  thickness  of  the  frozen  soil  layer. 

The  energy  extracted  (freezing)  or  available  (thawing)  at  a  node  ( J/m 3)  during  a  time  step 
can  be  expressed  as 

node  energy  =  A Qx  =  \Tt- 273. 15| ^Okpkcpt  (6. 1) 

where  7)  is  the  node  temperature  (K),  and  the  subscript  k  is  for  the  soil  components  dirt, 
water,  ice  and  air,  Ou  is  the  volume  fraction  of  component  k,  pk  ( kg/m  )  is  the  density  of 
component  k,  and  cpt  ( J/kg )  is  the  specific  heat  of  component  k.  Calculation  of  T,  is 

discussed  in  Chapter  4.  Determination  of  the  various  volume  fractions  is  discussed  in 
Chapter  5. 

The  energy  remaining  before  the  node  can  completely  freeze  or  thaw  (J/m3)  is  defined  as 
freeze  energy  =  A Q2  =  lfus  (0W  -  0r )  £  0kpkck  ^ 

thaw  energy  =  A  Q2  =  -/.  fus0^9kpkck 

where  lfus  is  the  latent  heat  of  fusion  of  water  (3.335  x  10 5  J/kg)  and  6,-  is  the  node  residu¬ 
al  volumetric  water  content  as  discussed  in  Chapter  7.  Freezing  can  occur  only  if  the  node 
f  <  273.15  K  and  thawing  if  Tt  >  273 . 1 5  if. 


The  change  in  volumetric  ice  content  per  time  step  is  therefore 
A0  _  Pw  min(A£„|Ag2|) 

Pi  lfus 

and  the  nodal  volumetric  soil  moisture  and  ice  contents  are  updated  accordingly. 
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6.2  Validation 


We  tested  the  freeze/thaw  module  using  the  data  collected  by  Lindamae  Peck  in  South 
Royalton,  VT,  spanning  the  winter  of  1990-1991  (personal  communication).  The  site  had 
silty  soil  of  which  the  density,  porosity,  and  thermal  conductivity  had  been  measured 
(Peck  and  O’Neill  1997).  Figure  8.1  compares  the  frost  depth  as  extrapolated  from  the 
thermister  probe  to  those  calculated  by  FASST.  FASST  was  initialized  with  two  soil 
layers,  the  top  one  a  silty  sand  (SM)  0.10  m  thick  and  the  bottom  layer  an  organic  silty 
clay  (OL)  1.40  m  thick.  The  FASST  determined  surface  temperature  denoted  “st  FASST” 
is  also  shown  in  Figure  8.1.  The  largest  differences  between  the  calculations  and  the 
measurements  occur  during  the  thawing  period.  FASST  partially  thaws  a  few  days  before 
the  measured  values.  This  is  in  part  due  to  differences  in  the  snow  meltout  date  between 
the  model  and  the  observations  as  discussed  below.  Shifting  of  the  surface  probe  during 
the  winter  is  also  to  blame. 


South  Royalton,  VT  1990-91 


| - measured - FASST - st  FASST  | 


0.9  n 


33178  33193  33208  33223  33238  33253  33268  33283  33298  33313  33328  33343  33358 

Julian  Day 


Figure  6.1  Comparison  of  frost  depth  between  FASST  and  measured  values  for  South 
Royalton,  VT. 


Figure  6.2  shows  that  several  daily  shallow  freeze/thaw  events  happen  before  and  after 
the  first  snow.  This  is  a  common  occurrence  in  temperate  climates.  From  Figure  6.2  one 
can  also  discern  that  when  there  is  snow  present  (measured  values),  FASST  accurately 
predicts  the  frost  depth.  Differences  between  FASST  and  measured  values  can  be  due  to 
the  thermister  string  spacing.  Another  reason  could  be  how  “frozen  soil”  is  defined 
between  the  two.  FASST  considers  a  node  frozen  if  there  is  any  ice  present  while  the 
measured  value  is  strictly  defined  by  temperature.  FASST  does  a  fairly  good  job  at 
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predicting  frost  depth,  especially  their  occurrence  and  under  a  snow  cover.  It  also  does  a 
good  job  at  reproducing  the  diurnal  freeze/thaw  cycles  that  occur  in  the  spring  and  fall. 
These  are  very  important  to  mobility  considerations  as  they  result  in  slippery  surface 
conditions. 

South  Royalton,  VT  1990-91 


| ^—measured  FASST  —snow] 


33175  33195  33215  33235  33255  33275  33295  33315  33335  33355 


Julian  Day 

Figure  6.2  Comparison  of  FASST  calculated  frost  depth  and  measured  values  as  a 
function  of  soil  depth  and  snow  depth. 
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Chapter  7 

Snow  Accretion,  Depletion,  and  Meltwater  Outflow  Model 


7.0  Introduction 

The  Snow  Accretion-Depletion  Module  uses  the  Snowmelt  Numerical  Analysis  Package 
(SNAP)  developed  by  Albert  and  Krajeski  (1998).  It  is  a  physically  based  approach  to 
modeling  snowmelt,  where  the  physics  of  flow  through  snow  are  considered  and  the  melt 
is  driven  by  an  energy  budget  at  the  snow  surface.  Snow  accretion  occurs  when  a 
snowfall  amount  is  given  in  the  input  or  when  precipitation  is  given  in  the  input  and  the 
input  air  temperature  is  below  freezing.  In  the  latter  case,  the  precipitation  amount  is 
converted  to  a  snowfall  amount.  Sections  7.1  and  7.2  are  from  Albert  and  Krajeski 
(1998). 


7.1  Governing  Equations  of  Flow  within  the  Snow 


For  modeling  the  movement  of  water  through  the  snow,  the  effects  of  capilarity  are  taken 
as  negligibly  small  compared  to  the  effects  of  gravity  (Colbeck  1972),  yielding  the 
simplified  form  of  Darcy’s  equation: 

PWK§ 


U  = 


(7.1) 


where  U  is  the  volume  flux  of  water,  pw  is  the  density  of  water,  kw  is  the  relative 
permeability  to  water,  g  is  the  acceleration  of  gravity,  and  q „  is  the  viscosity  of  water. 
Under  some  circumstances  this  will  be  applicable  to  the  entire  snowpack,  while 
modifications  will  be  necessary  under  some  conditions  of  layering.  The  effective 
permeability  of  the  water  phase  is  taken  to  be  proportional  to  a  power  ( nc )  of  the  effective 
water  saturation  (Morel-Saytoux  1969,  Colbeck  1972): 

K=K»sn;  (7.2) 


where  the  effective  water  saturation,  Se,  is  defined  by 


5.  = 


S  —  S 

uwa  wi 

1  -S', 


(7.3) 


where  S .  is  the  absolute  water  saturation  and  S.  is  the  irreducible  water  saturation.  The 

wa  wi 


general  applicability  of  the  relationship  kw  oc  S’’;  is  discussed  in  Maulem  (1978).  Here,  nc 
is  taken  as  a  constant  with  a  default  value  of  3.3. 


The  flux  in  terms  of  absolute  water  saturation  is  defined  as 


U  = 


Pw^woS 

vw 


(7.4) 


The  water  volume  conservation  equation  (Colbeck  1972)  used  states  that  the  change  in 
water  volume  flux  with  depth,  U  (where  the  z  axis  is  positive  downwards),  is  equal  to  the 
change  in  water  saturation  with  time: 
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U'  +  <f>(\-Swi)Se  =  0  (7.5) 

where  (f>  is  the  total  pore  volume  as  a  function  of  the  total  volume.  Changes  in  Sw,  occur 
on  a  much  larger  time  scale  than  changes  in  effective  saturation,  hence  it  is  sufficient  to 
consider  Swi  constant  over  an  individual  time  step  (although  it  can  be  altered  through  the 
course  of  the  snowmelt  season).  To  obtain  an  analytical  expression  that  can  easily  be 
employed  in  a  computer  model,  the  solution  of  these  equations  differs  from  earlier 
solutions  (Colbeck  1972,  Tucker  and  Colbeck  1977)  in  that  Equation  (7.4)  is  used  to  find 
Se  as  a  function  of  U,  which  is  then  substituted  into  Equation  (7.5)  above,  to  give  the 


governing  equation  for  water  volume  flux 


U  =  -ncf\l-Swiy 


1  f  Px\kwQ§ 


Uy 


In  order  to  make  the  problem  tractable,  nc,  cj),  Swi,  and  kwo  are  assumed  to  be  constant  over 
each  time  step.  These  variables  change  slowly  compared  to  the  time  scale  with  which 
mobile  water  moves  through  the  snow.  By  simplifying  notation  as  follows 


Kcf=ncfl{\-Swi)- 


-i  i  pwKoS 


the  flux  governing  equation  becomes 


U  =  -KcfU'lfU' . 


This  equation  is  solved  using  separation  of  variables.  The  general  solutions  to  the  volume 
flux  equation  are  therefore  sums  in  space  and  time  of  particular  solutions  governed  by 
boundary  conditions  while  maintaining  the  restrictions  indicated  above.  It  is  not 
necessary  to  completely  solve  the  volume  flux  equation,  as  a  temporal  expansion  of 
particular  solutions  will  tend  to  approximate  the  exact  solution  if  the  particular  solutions 
are  generated  near  enough  to  each  other  in  time.  With  each  separate  particular  solution 
being  generated  by  meteorological  data  in  separate  time  steps,  the  accuracy  of  the 
solution  will  be  dependent  on  the  time  step  of  the  available  meteorological  data,  with 
hourly  data  proving  accurate  enough  to  give  promising  results,  and  finer  resolution  data 
resulting  in  more  accurate  approximations  of  the  exact  solution. 


Equation  (7.9)  is  solved  by  assuming  the  solutions  take  the  fonn 
U(x,t)  =  f?  (x  +  q )a"  (t  +  c2Y 

where  substitution  back  into  Equation  (7.9)  and  re-grouping  gives 

p  (x + c2  r  (t+c2  r = -k^  acc  (x + q  r^-1  {t+C2  . 

For  this  to  hold 


a.„  =  — 


P  =  ~lt 


(7.10) 


(7.11) 


(7.12) 

(7.13) 


B  =  k  , 


(7.14) 
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(7.15) 


implying  the  solution  must  have  the  fonn 

U (x, t )  =  k^c{  (x  +  q )/C  [t  +  c2 )Af . 
Co-ordinates  are  then  nonnalized  by  defining 
X  =  x  +  q 


T  =  t  +  C2 

thus  giving 


(  X 


v  K<fT  j 


(7.16) 

(7.17) 

(7.18) 


In  the  nonnalized  co-ordinates  (%,  r),  each  particular  solution  represents  a  wave  of 
volume  flux  with  its  origin  at  the  point  ( x  =  0,  r  =  0  ).  In  other  words,  each  differential 

element  of  mobile  water,  dAm,  generates  its  own  particular  solution  where  t  =  0  occurs  at 
the  time  the  water  volume  becomes  mobile  within  the  snow  (surface  melt  or  rainfall),  and 
X  =  0  occurs  at  the  vertical  height  of  the  snowpack  at  the  time  the  water  volume  becomes 
mobile. 


To  find  the  depth  of  penetration  of  a  wave  as  a  function  of  time  is  simply  a  matter  of 
integrating  saturation  from  the  point  of  generation  of  the  wave  to  whatever  depth 
necessary  to  allow  the  total  volume  enclosed  by  the  integration  to  be  equal  to  the  original 
volume  of  the  wave.  For  an  input  volume  per  unit  area  of  A  (cm)  over  an  area  <jwa,  if  the 
depth  of  penetration  of  the  resulting  flux  wave  is  denoted  as  xu,  then  at  any  point  t  in 
time 


f  XJ  S/M]  ~  Swi)(Jwadx  = 

•'7=0 


‘'7=0 

where  A  is  the  volume/unit  area  runoff  in  cm  and  <jwa  is  the  area.  Therefore  we  have 


(7.19) 


J 


"¥z=— 

*-°  «Ki-s„) 


(7.20) 


Using  Se  from  Equation  (7.4)  and  U  from  Equation  (7.18),  substituting  into  Equation 
(7.20),  then  solving  the  integral  yields  the  equation  for  the  depth  of  penetration  of  a  wave 
as  a  function  of  time 

r  ■  vrDk 

PwKw  OS 


xv  = 


A 


Md  -sj 


V 


) 


(*VT) 


(7.21) 


\  I  cjr  ~wi  /  j 

Next,  consider  the  time  it  will  take  a  wave  of  volume  flux  to  penetrate  to  a  given  depth 
(most  useful  in  calculating  the  delay  between  influx  and  outflux).  In  order  for  the  water 
mass  to  balance  in  time,  the  volume  of  water  passing  by  a  given  depth  D  must  equal  the 
initial  mobile  water  volume 


f  U(D,r)dr  =  Am 

JT=Tn 


Substituting  from  Equation  (7.15)  yields 

-i  i  -i 

r^D^Pdr  =  A„ 


where  D  is  the  depth  of  penetration  at  time  To- 


(7.22) 

(7.23) 
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The  solution  to  this  integral  gives  the  time  required  for  a  wave  of  volume  flux  to 
penetrate  to  depth  D 


^0  = 


l~Vcf 

'If 


A 


K_4_ 

\D  J 


% 


if 

Vcf-1 


(7.24) 


) 


For  modeling  purposes,  note  that  Equation  (7.21)  may  also  be  written  in  the  simpler  form 

(  n^-l 
v  A 


XU  =  Kcf 


iVr-  - 


r  c  . 


(7.25) 


J 


Having  depth  of  penetration  as  a  function  of  time,  and  vice  versa,  requires  only  that  the 
effects  of  preceding  flux  waves  be  incorporated  in  order  to  construct  a  working  model  of 
water  movement  through  homogeneous  snow.  Therefore,  consider  two  flux  waves,  each 
proceeding  according  to  Equation  (7.18) 

i 


Z_ 

KKTXJ 


If 


(7.26) 


and 


^2  {Z2’Tl) 


Z2 

V  KcfZ 2  J 


If 


(7.27) 


Assume  that  Xi  -  =  £  >  0  (that  the  wave  U\  began  at  an  earlier  point  in  time  than  the 

wave  U2).  Since  there  is  no  singularity  within  the  region  over  which  the  U2  solution 
applies,  flux  will  be  continuous  and  differentiable  over  this  region,  up  to  the  singularity 
caused  by  the  U\  -  Exjunction.  As  the  U2  wave  flows  through  the  snow,  Equation  (7.25) 
must  always  be  met,  even  as  new  water  volume  is  absorbed  into  wave  U2  from  the 
residual  saturation  of  wave  U\.  Therefore,  to  accommodate  the  new  water  volume  the 
second  wave  has  encountered  while  still  maintaining  its  solution  to  the  conservation  and 
flux-saturation  equations,  the  volume  flux  wave  EX  must  travel  deeper  into  the  snowpack 
than  it  otherwise  would,  thereby  encountering  more  residual  saturation  in  the  process. 
Numerically,  the  residual  saturation  encountered  will  take  the  form 


S,= 


>1, 

Pw^woS 


1 

nc 


(7.28) 


which,  when  combined  with  U\  from  Equation  (7.26),  becomes 

1  1 


7, 

nc 

(  Zl 

_Pw^woS  \ 

l*VrJ 

(7.29) 


which  means  that  if  xai  is  the  deepest  point  within  the  snow  at  which  the  EX  solution 
applies,  (the  point  of  the  singularity),  then  the  volume  of  water  absorbed  by  the  EX  wave 
in  moving  some  infinitesimal  distance  through  the  snow  will  be  given  by 

1 


7, 

nc 

(  \ 
XA2 

_Pw^wO§  _ 

-1 


dx 


A  2  ' 


(7.30) 
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Therefore,  as  found  before,  (derived  with  constant  volume,  but  still  applicable  with  a  time 
varying  volume), 


1 — « 


VA2 


fi)  = 


K 


cf 


-1 


V  A2  J 


(7.31) 


where  now  A2  is  a  function  of  time.  Therefore,  define  a  constant  An  equal  to  the  initial 
mobile  water  volume  of  the  U2  wave,  and  a  new  variable  A?  equal  to  the  total  water 
volume  the  U2  wave  has  absorbed,  by  the  equation 

A2  =  Au  +  A2  (7.32) 

where  we  know  that  because  the  U\  wave  is  undisturbed  beneath  the  U\  -  f/2  junction,  the 
total  water  volume  below  the  junction  must  be  the  same  as  it  would  be  were  the  U2  wave 
not  present,  giving  rise  to  the  conclusion  that  the  total  water  volume  above  the  U1-U2 
junction  must  be  equal  to  the  initial  volume  of  the  U2  wave  plus  the  volume  of  the  U\ 
wave,  which  would  not  otherwise  have  passed  the  U1—U2  junction  point.  This  means  that 


A2  = 


A2 


K 


cf 


(«c-iR"c- 


By  substituting  into  Equation  (7.32)  we  get 


(7.33) 


Xa2{t2)  = 


K 


cf 


n  - 1 


An  + 


XA2 
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cf 


(«c -l)(^i)‘"c 


r2‘ , 


which  can  be  solved  forx^2 


XA  2  Kcf^2 


(  A  7 
Au 

nc- 1 

1- 

1 

nc- 1 

lnc~lJ 

(7.34) 


(7.35) 


This  is  the  fonn  used  in  the  model  to  find  the  depth  of  any  wave  that  has  a  preceding 
wave.  Finding  the  point  at  which  one  wave  of  volume  flux  completely  overtakes  its 
antecedent  proves  to  be  straightforward  if  the  first  wave  has  no  predecessor,  with  the 
result  being 

(7.36) 


r2=S 


f  a 

^IL+\ 


\  Ai 


-1 
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where  8  is  the  time  between  the  fonnation  of  wave  U\  and  wave  U2.  This  gives  us  the 
ability  to  calculate  the  flux  as  a  function  of  depth  and  time  resulting  from  a  single  water 
volume  input,  to  calculate  the  depth  and  net  water  volume  of  a  following  input  wave  of 
volume  flux,  and  to  calculate  the  depth  and  time  at  which  one  wave  of  volume  flux  will 
overrun  a  preceding  wave  of  volume  flux.  In  short,  it  is  now  possible  to  completely 
describe  the  saturation  and  flux  resulting  from  two  point  inputs  in  several  compact 
equations,  without  any  finite  difference  or  finite  element  matrix  solutions  to  differential 
equations.  In  order  to  extend  the  model  to  the  general  case  in  which  many  waves  may  be 
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flowing  through  the  snowpack  without  allowing  time  for  complete  drainage  or  refreezing 
between  pairs  of  volume  flux  waves,  however,  it  is  necessary  to  know  the  time  of 
collision  of  a  wave  of  volume  flux  with  its  predecessor,  which  itself  is  not  traveling 
through  new  snow.  Unfortunately,  in  the  event  that  there  are  multiple  waves  flowing 
through  the  snow,  no  explicit  form  for  the  time  of  collision  can  be  found.  Furthermore,  no 
explicit  equation  for  the  time  a  given  volume  flux  will  take  to  reach  any  given  depth  can 
be  found.  This  is  not  the  downfall  of  the  model,  however,  because,  while  the  equations 
governing  the  time  of  collision  between  waves  are  not  tractable  in  general,  they  are  very 
well  behaved,  becoming  linear  very  shortly  after  the  generation  of  the  wave.  Therefore 
the  snow  melt  module  uses  a  Newton’s  method  approximation  to  find  the  time  a  given 
volume  flux  will  take  to  reach  the  bottom  of  the  snowpack  (and,  if  necessary,  the  time  of 
collision  of  two  waves  in  the  general  case).  This  usually  takes  four  or  less  iterations  of 
the  approximating  loop  to  come  within  ±  0.5  %  of  the  true  value. 

The  process  of  water  volume  flux  waves  absorbing  residual  mobile  water  saturation  from 
the  tails  of  their  predecessors  (described  by  Equation  [7.36]) ,  which  states  that  any  wave 
will  eventually  overtake  any  preceding  wave  given  an  infinite  time  and  an  infinite  snow 
depth  to  traverse),  although  it  necessitates  the  Newton’s  method  approximation  above, 
also  leads  to  the  great  strength  of  this  model:  its  computational  simplicity.  Because  waves 
continually  overtake  each  other  on  the  way  to  the  bottom  of  the  pack,  there  rarely  are 
very  many  of  them  within  the  snowpack  at  any  given  time.  The  only  conditions  under 
which  there  will  be  very  many  waves  within  the  snowpack  is  when  dAm  is  low  enough, 
and  snow  depth  is  high  enough  that  melt  waves  take  a  long  time  (compared  to  the  time 
scale)  to  flow  out  of  the  pack,  and  input  volume  flux  is  decreasing  at  a  rapid  enough  rate 
that  waves  are  not  colliding  within  the  snow.  Using  an  hourly  time  step,  and  a  surface 
energy  balance  for  melt  input,  it  is  rare  to  see  more  than  three  waves  within  the  snow  at 
any  one  time  in  the  seasonal  snow  of  New  England.  Any  further  modeling  of  flux  waves 
on  their  way  through  the  pack  would  not  result  in  any  greater  accuracy  of  outflow 
prediction,  as  only  one  wave  may  actually  flow  out  of  the  bottom  of  the  pack  at  any  time. 


7.2  Evolution  of  Parameters 

While  the  snow  property  parameters,  such  as  grain  size,  kwo,  Swi,  nc,  and  (jy  are  relatively 
stable  over  the  usual  lifetime  of  a  wave  of  volume  flux,  over  the  lifetime  of  a  snowcover 
they  do  vary.  Here  the  parameters  are  considered  as  constants  in  a  given  time  step,  yet 
they  may  change  over  the  course  of  the  snowmelt  season. 


7.2.1  Grain  Growth 

Grain  size  is  currently  used  only  to  calculate  the  permeability  of  snow.  The  equations  of 
Brun  (1989)  are  used  to  calculate  average  grain  volume.  The  initial  volume  of  the  snow 
crystals  is  currently  taken  as  0.00005  cm3,  which  is  quickly  dwarfed  by  the  growth 
equations,  especially  if  the  snow  is  wet.  Once  the  average  crystal  volume  in  the  wet 
portion  of  the  pack  and  the  average  crystal  volume  in  the  dry  portion  of  the  pack  are 
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calculated,  a  weighted  average  of  the  two  is  taken  as  the  average  crystal  volume  within 
the  pack. 


7.2.2  Permeability 

The  most  exhaustive  search  published  to  date  for  the  permeability  of  snow  as  a  function 
of  grain  size  and  snow  density  probably  comes  from  the  work  of  Shimizu  (1970),  in 
which  he  considered  previous  models  of  snow  penneability  produced  by  Bader  (1954), 
Bender  (1957),  and  Kuroiwa  (1968),  as  well  as  a  range  of  his  own  tests,  and  came  to  the 
conclusion 

kw0  =  0.077 d2  exp(-7.8/?s)  (7.37) 

where  d  is  the  snow  grain  diameter  (cm.)  and  pv  is  the  specific  gravity  ( g/cm ).  In  our 
model,  d  is  computed  as 

<7-38> 

where  Vav  is  the  average  crystal  volume  computed  above.  p.s  is  taken  as  the  density  of 
frozen  water  within  the  snowpack,  either  as  computed  in  the  snow  depth  prediction 
routine,  or  as  a  constant  0.3  g/cm  if  snow  depth  is  not  being  computed. 

7.2.3  Irreducible  Water  Saturation 

In  the  FASST  snow  module,  the  irreducible  water  saturation  of  snow  is  computed  as  a 
fraction  of  total  volume.  Kattlemann  (1986)  reported  that  Sulahrie  (1972)  found  a  value 
as  high  as  40%  in  some  cases,  while  most  authors  have  reported  measurements  in  the 
range  of  0%  to  10%. 


7.2.4  Snow  Depth 


Snow  depth  prediction  equations  are  based  on  the  form  used  by  Jordan  (1991a),  which  in 
turn  come  from  Anderson  (1973).  Both  SNTHERM.89  (Jordan  1991a)  and  Anderson’s 
earlier  model  break  the  snow  into  layers  based  on  age,  density,  and  crystal  structure, 
while  here  the  simplicity  of  a  uniform  bulk  snowpack  is  maintained.  The  model  predicts 
the  rate  of  densification,  then  adjusts  the  snowdepth  accordingly.  The  equations  used  are 
as  follows 


1  PD. 


D.  dt 


=  -2.778  x \W6cxc2  exp[-0.04T] 


(7.39) 


metamorphism 

where  Ds  is  the  depth  of  snow,  t  is  time  (5),  T  is  temperature  (°C),  and 

Cj  =  1  pt  <  0.12 


q  =  exp  -46 (ft- 0. 1 5)] 


Pi  -  0-12 


(7.40) 


where  p,  is  the  density  of  water  in  the  frozen  state  within  the  snowpack  (g/cm.)  and 
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(7.41) 


c2  =  1  +  // 

where  f  is  the  fraction  of  the  snowpack  which  is  wet.  Further,  we  have 


1  dD. 


D,  dt 


A 

*ic 


(7.42) 


\  overburden 

where  Ps  is  the  average  load  pressure  within  the  snowpack,  and  rjc  is  a  viscosity 
coefficient,  which  has  the  form 

7e  =7o  exp (0.087  +  28/?,)  (7.43) 

where  p,  is  the  total  density  of  the  solid  and  liquid  phases  of  the  snow,  and  r|o  =  5E+08. 


Also  calculated  with  snow  depth  are  /?,■  (the  density  of  ice  within  the  snowpack),  p,  (the 
density  of  ice  plus  irreducible  water  saturation  plus  mobile  water  saturation),  <j>  (the  total 
pore  volume  as  a  fraction  of  total  volume),  (j)e  (pore  volume  excluding  pore  volume  filled 
with  immobile  water  when  meeting  irreducible  saturation  requirements,  as  a  fraction  of 
total  volume),  and  SWE  (the  snow  water  equivalent). 


7.2.5  Effective  Saturation  Exponent 

The  general  applicability  of  the  relationship  kw  oc  S"c  is  discussed  in  depth  in  Mualem 
(1978).  Since  a  way  to  determine  nc  from  the  available  meteorological  and  lysimeter  data 
could  not  be  found,  FASST  simply  uses  a  constant  nc,  and  its  value  is  left  up  to  the  user, 
(with  a  default  of  3.3). 


7.2.6  Refreezing 

The  refreezing  algorithm  uses  a  time-averaged  value  of  temperature  over  the  most  recent 
period  in  which  the  snowpack  is  predicted  to  be  less  than  isothermal,  and  a  depth- 
averaged  value  of  saturation,  along  with  bulk-approximated  thermal  conductivity  of 
snow,  to  calculate  the  depth  of  penetration  of  the  refreezing  front.  Since  these  values  are 
updated  every  time  step,  the  depth  of  penetration  produced  should  be  accurate  enough. 
An  analytical  solution  of  the  well-known  Neumann  type  (Carslaw  and  Jaeger  1959)  given 
by 


X  = 


2k\Tf-T)t 


Plf 


(7.44) 


where  X  is  the  depth  of  the  freezing  front,  k\  is  the  thermal  conductivity  of  the  medium  in 
the  refrozen  state,  Tf  is  the  temperature  of  fusion,  Ts  is  the  surface  (ambient  air) 
temperature,  t  is  time,  p  is  the  density  of  the  medium,  and  //is  the  latent  heat  of  fusion. 


k\  is  taken  to  be  0.0045  J/s  cm-°C,  which  is  approximately  the  thermal  conductivity  of 
snow  of  density  0.3  g/cm3.  Tf  is  taken  as  0  °C,  while  Ts  is  taken  as  the  average 
temperature  over  the  current  period  of  freeze  depth  propagation,  //is  taken  as  333.05  J/g, 
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so  that  the  latent  heat  of  fusion  of  1  cm 3  of  snow  that  has  been  drained  to  Swi  will  be  l-Swi. 
A  correction  to  account  for  residual  mobile  water  saturation  was  implemented,  but  the 
effects  were  found  to  be  minimal  in  light  of  the  fact  that  the  freezing  front  progresses 
slowly  compared  to  the  speed  of  water  percolation. 


7.3  Surface  Energy  Balance 


In  generating  melt,  FASST  currently  uses  a  full  surface  energy  balance  to  calculate  the 
volume  of  runoff  that  is  generated  from  surface  melt  during  each  time  step.  In  general, 
the  heat  input  at  the  top  of  the  snow  (I top)  is  (Jordan  1991a) 

l,v=lJ(i-a)  +  li-ll  +  H  +  L  +  (7.45) 

where  Is  X  is  the  net  solar  radiation  at  the  surface,  a  is  the  surface  albedo  pertinent  to 

I  * 

shortwave  radiation,  Iir  is  the  incoming  longwave  radiation,  Iir  is  the  outgoing  longwave 

radiation,  H  is  the  sensible  heat  flux,  L  is  the  latent  heat  flux,  and  is  the  convective 
heat  flux.  Solar  radiation  and  albedo  are  currently  not  estimated  in  the  module,  and 

neither  is  net  longwave  radiation,  so  Is  i  ( 1  -  a)  and  [lfr  -  ljr  j  are  taken  as  inputs  to  the 

model.  Details  of  the  exchange  follow  Jordan  (1991a,  b).  The  resulting  surface  melt  depth 
is 

At 


Ahi,<oP  =  hop 


Pih 


(7.46) 


If  the  ground  temperature  is  above  freezing,  melting  can  also  occur  at  the  bottom  of  the 
snow.  The  bottom  melt  depth  is 

a  /  dT  At  . 

Ahbottom  =  - Y  (7-47) 

dz  p  i  , 

where  k  is  the  soil  thermal  conductivity  ( W/m-K ),  T  is  the  soil  temperature  ( K ),  and  z  is 
the  soil  depth  (m). 


7.4  Module  Input  and  Output 

The  Snow  Accretion/Depletion  Module  requires  certain  information  about  the  snow 
cover.  If  known,  snow  depth,  snow  water  equivalent,  initial  water  saturation  (default 
0.0288),  effective  porosity,  and  the  snow  surface  temperature  should  be  input.  Otherwise, 
the  module  will  compute  these  parameters.  The  required  meteorological  data  are  air 
temperature,  wind  speed,  net  solar  radiation,  net  longwave  radiation,  and  precipitation 
(amount  and  type).  If  the  type  of  precipitation  (rain,  snow)  is  not  known,  it  will  be 
predicted  based  on  the  ambient  air  temperature.  Net  radiation  can  be  obtained  using  the 
two  radiation  modules  described  in  Chapters  4  and  5.  The  output  of  the  module  is  snow 
depth  and  the  volume  of  water  (total  and  incremental)  that  has  come  out  of  the  snowpack. 
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7.5  Model  Validation 


To  test  how  well  FASST  predicts  the  snow  depth,  we  compared  the  model  results  with 
data  gathered  during  the  winter  of  1997  in  the  Sleepers  River,  VT,  watershed  and  during 
the  winter  of  2003  for  Buffalo  Pass,  CO,  and  Indian  River,  CO.  The  latter  two  data  sets 
are  part  of  NASA’s  Cold  Land  Processes  Experiment  (CLPX).  The  results  for  Sleepers 
River  are  plotted  in  Figure  7.1.  FASST  can  run  in  two  modes,  one  where  FASST  uses  the 
measured  snow  depth  to  self-correct  the  model  predictions  (“FASST,  Yes”)  and  one 
where  only  initial  snow  depth  is  provided  (“FASST,  No”).  The  relative  errors  are  1%  and 
7%,  respectively,  both  well  within  the  accepted  accuracy  standard  of  15%  (Flolcombe  et 
al.  2004).  It  can  also  be  seen  in  Figure  7.1  that,  for  the  most  part,  increases  in  snow  depth 
are  accompanied  by  precipitation  events.  Increases  in  the  measured  values  not  associated 
with  snowfall  are  the  result  of  blowing  snow.  FASST  does  not  account  for  this  at  present. 


Sleepers  River,  VT  1997 


80  85  90  95  100  105  110  115  120 

Day  of  the  Year 


Figure  7.1  Comparison  between  the  measured  snow  depth  and  the  snow  depth  calculated 
by  FASST  along  with  snow  (yellow)  and  rain  (blue)  precipitation  events  for  Sleepers 
River,  VT. 
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Buffalo  Pass 


| - measured - sntherm - fasst] 


Figure  7.2  Comparison  between  the  measured  snow  depth  and  the  snow  depth  calculated 
by  SNTHERM  and  FASST  for  Buffalo  Pass,  CO. 

Illinois  River 


50  60  70  80  90  100  110  120  130 


day  of  the  year 

Figure  7.3  Comparison  between  the  measured  snow  depth  and  the  snow  depth  calculated 
by  SNTHERM  and  FASST  for  Indian  River,  CO. 
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Figures  7.2  and  7.3  compare  the  results  of  FASST  with  those  of  SNTHERM  against  the 
measured  values  for  two  sites  in  Colorado.  SNTHERM  (Jordan  1991a)  is  considered  the 
most  physically  based  snow  model.  The  Buffalo  Pass  site  is  characterized  by  an 
accumulation  period  followed  by  an  ablation  period  that  begins  on  approximately  day 
130.  Both  models  accurately  transition  between  the  two  seasons.  FASST  underestimates 
the  meltout  date  by  three  days  while  SNTHERM  overestimates  it  by  three  days.  Within 
five  days  is  considered  accurate  (Holcombe  et  al.  2004).  The  absolute  (average  of  model- 
measured)  and  relative  (average  of  [model-measured]/measured)  errors  for  both  models 
are  8%  and  3%,  respectively. 

The  snowpack  at  the  Illinois  River  site,  Figure  7.3,  is  characterized  as  thin,  ephemeral, 
and  windblown.  This  is  a  more  difficult  modeling  scenario.  Despite  this,  both  models  do 
well.  The  absolute  error  is  <  1%  while  the  relative  error  is  2  %  for  SNTHERM  and  11% 
for  FASST. 
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Chapter  8 

Ice  Thickness  Model 


8.0  Introduction 

Of  great  hindrance  to  mobility  is  the  presence  of  ice  on  the  ground.  This  can  happen 
when  the  ground  is  at  or  below  freezing  and  any  precipitation  is  in  the  form  of  rain 
instead  of  snow. 


8.1  Theory 

8.1.1  Growth 


The  formation  of  an  ice  cover  on  the  surface  is  allowed  only  if  there  is  no  snow.  If  snow 
is  present,  any  rain  that  falls  is  assumed  to  percolate  into  the  snowpack  instead  of  forming 
an  ice  cover  on  the  snow  surface.  Following  Jones  (1996),  the  fraction  of  precipitation 
that  freezes  is 


Jheat  fluxes -latent  heat  flux  (8  1) 

latent  heat  flux 

The  sum  of  the  heat  fluxes  is  fully  described  in  Chapter  6,  Section  6.1.  The  resulting  ice 
thickness,  ht  (in)  per  unit  area  is 


h,  =  fP - ^  (8.2) 

'  3600  Pi 

where  P  (m/hr)  is  the  precipitation  rate,  At  (sec)  is  the  time  step,  pw  (kg/m3)  is  the  density 
of  water,  and  p,  (kg/m)  is  the  ice  density. 


8.1.2  Decay 


Decay  of  an  ice  layer  can  occur  from  both  the  top  and  bottom,  but  only  if  enough  energy 
is  present  at  the  corresponding  interface.  Once  this  condition  is  met,  the  resulting  melt 
depths,  Ahi'top  (m)  and  Ap, bottom  (m),  are 


Ahi,op  =  Vnet 


AhLbo,tom  =  K 


At 
Pjf 
dt  At 
dz  ptlf 


(8.3) 


2 

where  qnet  (W/m  )  is  the  net  surface  heat  flux,  f  (J/kg)  is  the  latent  heat  of  fusion, 
Ci  (J/kg-K)  is  the  specific  heat  of  ice,  Ts  (K)  is  the  surface  temperature,  and  k  is  the 
soil  thermal  conductivity  (W/m-K).  This  is  the  same  procedure  used  by  the  snow 
accretion  and  depletion  model  described  in  Chapter  9. 
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Chapter  9 

Model  Inputs 


9.0  Introduction 

FASST  requires  four  types  of  inputs:  meteorological,  soil  properties,  initial  conditions, 
and  site  information.  Actual  structure  of  the  input  files  is  not  discussed  here.  This 
information  is  found  in  the  accompanying  report,  FASST  User’s  Manual  (Frankenstein 
2004).  Wherever  possible,  SEDRIS  (Synthetic  Environment  Data  Representation  and 
Interchange  Specification)  nomenclature  is  used  to  describe  the  data  in  an  attempt  to 
standardize  the  process.  More  information  concerning  SEDRIS  is  found  at 
http://www.sedris.org. 


9.1  Meteorological  Data 

For  FASST  to  operate  with  the  most  accuracy,  the  following  meteorological  data  are 
required. 

Year 

Day  of  the  Year 
Hour  of  the  Day 
Minute  of  the  Hour 
Air  Pressure  ( mbar ) 

Air  Temperature  ( K) 

Relative  Humidity  (%) 

Wind  Speed  (m/s) 

Wind  Direction  ( degrees  from  N,  +  =  clockwise) 

Precipitation  Rate  (mm/hr) 

Precipitation  Type  (SEDRIS  enumeration) 

Low  Cloud  Amount  ( tenths ) 

Low  Cloud  Height  (km) 

Low  Cloud  Type  (SEDRIS  enumeration) 

Middle  Cloud  Amount  (tenths) 

Middle  Cloud  Height  (km) 

Middle  Cloud  Type  (SEDRIS  enumeration) 

High  Cloud  Amount  (tenths) 

High  Cloud  Height  (km) 

High  Cloud  Type  (SEDRIS  enumeration) 

Total  Incoming  Solar  Radiation  (Direct  +  Diffuse)  (W/m2) 

Direct  Solar  Radiation  (W/m2) 

Diffuse  Solar  Radiation  (W/m2) 

Reflected  Solar  Radiation  (W/m2) 

Net  Solar  Radiation  (Total  -Reflected)  (J/m2) 

Incoming  Infrared  Radiation  (W/m2) 
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2. 

Emitted  Infrared  Radiation  ( W/m  ) 

Net  Infrared  Radiation  ( J/m 2) 

Solar  Zenith  Angle  ( degrees ) 

Solar  Azimuth  Angle  ( degrees ) 

Snow  Depth  (m) 

Soil  Surface  Temperature  ( K) 

At  a  minimum,  FASST  requires  that  the  user  provide  the  observed  air  temperature  along 
with  the  time  and  date.  The  data  frequency  may  range  from  1  minute  to  24  hours,  with 
the  optimum  range  being  1-3  hours.  If  any  of  the  above  parameters  are  not  provided,  the 
following  assumptions  are  made. 


Air  Pressure 

1000  mbar 

Relative  Humidity 

60% 

Wind  Speed 

1 .2  m/s 

Wind  Direction 

0.0  degrees 

Precipitation  Rate 

0.0  mm/hr 

Precipitation  Type 

1  (none) 

Low  Cloud  Amount 

0.5  tenths 

Low  Cloud  Height 

1.5  km 

Low  Cloud  Type 

6  (stratus  nebulosus  or  stratus  fractus) 

Middle  Cloud  Amount 

0.0  tenths 

Middle  Cloud  Height 

4.0  km 

Middle  Cloud  Type 

3  altocumulus  translucidus,  1  level 

High  Cloud  Amount 

0.0  tenths 

High  Cloud  Height 

8.0  km 

High  Cloud  Type 

5  (cirrus  and/or  cirrostratus  <  45°  above  the  horizon) 
if  High  Cloud  Amount  <  0.4 

7  (cirrostratus,  full  cover)  if  High  Cloud  Amount  >  0.4 

Total  Incoming  Solar  Radiation  (Direct  +  Diffuse) 

Direct  Solar  Radiation 

Diffuse  Solar  Radiation 

Upwelling  Solar  Radiation 

Net  Solar  Radiation  (Total  -  Upwelling) 

Incoming  Infrared  Radiation 
Emitted  Infrared  Radiation 
Net  Infrared  Radiation 
Solar  Zenith  Angle 
Solar  Azimuth  Angle 
Snow  Depth 

Soil  Surface  Temperature 


Model  will  calculate 
Model  will  calculate 
Model  will  calculate 
Model  will  calculate 
Model  will  calculate 
Model  will  calculate 
Model  will  calculate 
Model  will  calculate 
Model  will  calculate 
Model  will  calculate 
Model  will  calculate 
Model  will  calculate 


The  enumerations  used  for  the  precipitation  and  cloud  types  are  based  on  SEDRIS. 
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9.2  Soil  Properties 

Besides  the  meteorological  data,  how  representative  the  soil  parameters  are  of  the  region 
being  modeled  has  a  significant  effect  on  FAS  ST ’s  ability  to  accurately  predict  the  state 
of  the  ground.  At  a  minimum,  the  user  must  provide  infonnation  concerning  the  number 
of  layers  (maximum  is  10),  layer  thickness,  and  the  layer  USCS  (Unified  Soil 
Classification  System)  soil  type.  The  soil  parameters  that  FASST  uses  are 


Bulk  density  of  dry  material  ( g/cm 3)  yd 

Intrinsic  density  of  dry  material  (g/cm3)  pd 

Volume  fraction  of  solids  Od 

Porosity  n 

Void  Ratio  e 

Albedo  (0.35-3.0 pm)  a 

Emissivity  (3.0-100  pm)  s 

Quartz  content  q 

Organic  Fraction  0oj 

Thermal  conductivity  of  dry  material  (W/m  ■ K)  kth 

Specific  heat  of  the  dry  material  ( J/kg-K)  C 

Saturated  hydraulic  conductivity  (cm/s)  Ksat 

Residual  water  content  (vol/vol)  0, 

Maximum  water  content  (vol/vol)  0max 

van  Genuchten  bubbling  pressure  head  (cm)  ccvg 

van  Genuchten  exponent  nVQ 

Cone  index/moisture  content  coefficient  1  a 

Cone  index/moisture  content  coefficient  2  c? 

Rating  cone  index/moisture  content  coefficient  1  a 
Rating  cone  index/moisture  content  coefficient  2  c? 


Many  of  the  above  parameters  are  related.  For  instance,  the  bulk  density  (ya),  mass  per 
total  volume,  is  related  to  intrinsic  density  (pd),  mass  per  fractional  volume  of  solids,  in 
the  following  manner  yd  =  9dpd  where  Od  is  the  volume  fraction  of  the  solids  to  the  soil 
as  a  whole.  Remember  also  that  the  porosity,  n,  and  void  ratio,  e,  are  related  through 
n=e/(\  +  ey  In  addition,  the  porosity  is  also  a  function  of  the  volume  of  solids,  Od,  i.e, 

n  =  1  -  Od. 

The  thermal  conductivity  of  the  bulk  dry  material  (tQry)  is  related  to  the  thermal 
conductivity  of  the  solids  ( tcs)  by 

(9.1) 

where  rca  =  0.026  W/m-K  is  the  thermal  conductivity  of  air  and  n  is  the  porosity. 
Unfortunately,  the  thennal  conductivity  of  the  solids  is  seldom  known.  A  different 
approach,  and  the  one  used  by  FASST,  is  to  calculate  Kdry  based  on  (Farouki  1981): 
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Kdry 


0.13Sr,+64.7 

2700-0.947/^ 


(9.2) 


The  thermal  conductivity  is  related  to  the  specific  heat  such  that  t<  =  kth  C  where  kth  is  the 
thermal  diffusivity  (nr/s)  and  C  is  the  specific  heat  ( J/kg-K ). 


Finally,  it  should  be  noted  that  often,  the  saturated  hydraulic  conductivity,  Ksat  (cm/s),  is 

2 

called  the  saturated  permeability,  kpsat  (cm  ).  The  two  are  related  by 

(9.3) 

v  rtjpf 

where  A  is  a  pore  shape  dependent  proportionality  factor,  d  (in)  is  the  soil  grain  diameter, 
g  (m/s  )  is  the  acceleration  due  to  gravity,  v  (nr/s)  is  the  kinematic  fluid  viscosity,  rj 
(kg/m-s)  is  the  dynamic  fluid  viscosity,  and  pf  is  the  fluid  density  (kg/m). 

FASST  provides  default  values  for  all  of  the  above  parameters,  for  15  of  the  USCS  soil 
types  (see  Tables  9.2.1  and  9.2.2).  Parameters  not  found  in  the  tables  are  calculated  based 
on  the  preceding  discussion.  The  USCS  soil  types  that  FASST  allows  are 


GW  Well-graded  gravels,  gravel-sand  mixtures,  little  or  no  fines 

GP  Poorly  graded  gravels,  gravel-sand  mixtures,  little  or  no  fines 

GM  Silty  gravels,  gravel-sand-silt  mixtures 

GC  Clayey  gravels,  gravel-sand-clay  mixtures 

SW  Well-graded  sands,  gravelly  sands,  little  or  no  fines 

SP  Poorly  graded  sands,  gravelly  sands,  little  or  no  fines 

SM  Silty  sands,  sand-silt  mixtures 

SC  Clayey  sands,  sand-clay  mixtures 

ML  Inorganic  silts  and  very  fine  sands,  rock  flour,  silty  or  clayey  fine 

sands  or  clayey  silts  with  slight  plasticity 
CL  Inorganic  clays  of  low  to  medium  plasticity,  gravelly  clays,  sandy 

clays,  silty  clays,  lean  clays 

OL  Organic  silts  and  organic  silty  clays  of  low  plasticity 

CH  Inorganic  silts,  micaceous  or  diatomaceous  fine  sandy  or  silty 

soils,  elastic  silts 

MH  Inorganic  clays  of  high  plasticity,  fat  clays 

OH  Organic  clays  of  medium  to  high  plasticity,  organic  silts 

Pt  Peat  and  other  highly  organic  soils 

SMSC  Mixed  SM  and  SC  soils 

CLML  Mixed  CL  and  ML  soils 

EV  Evaporites  (salt  pans,  similar  properties  to  SW) 

CO  Concrete 

AS  Asphalt 

RO  Bed  rock  (assumed  to  be  granitic) 

SN  Pennanent  snow  field  or  glacier 

UK  Unknown  (This  is  defaulted  to  SM.) 

US  User-supplied  variables 
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If  the  soil  type  “US”  is  used,  the  user  is  then  directed,  in  a  step-by-step  manner,  to  enter 
their  values  for  the  soil  parameters.  If  any  of  the  parameters  are  unknown,  FASST  reverts 
to  the  default  values  listed  in  Tables  9.2.1  and  9.2.2  as  stated  above. 


Table  9.2.1  FASST  Default  Soil  Properties 


Soil 

Type 

Bulk 

Dry 

Density 

Yd  (g/cm3) 

Porosity 

n 

Emissivity 1 
£ 

Albedo1 

a 

Quartz 
Content 2 

q 

Organic 

Fraction 

Oof 

Fine/ 

Coarse 

Specific 

Heat3 

C,  (J/kg-K) 

GW 

1.95s 

0.296s 

0.92 

0.40 

0.65 

0.00 

1 

820.0 

GP 

2.16s 

0.203s 

0.92 

0.40 

0.65 

0.00 

1 

820.0 

GM 

1.91s 

0.324s 

0.95 

0.40 

0.65 

0.00 

1 

820.0 

GC 

1.87s 

0.344 5 

0.92 

0.40 

0.65 

0.00 

1 

820.0 

SW 

■IKttliM 

0.92 

0.40 

0.80 

0.00 

1 

830.0 

SP 

1.594s 

0.415s 

0.92 

0.35 

0.80 

0.00 

1 

816.47 

SM 

1.474s 

0.526s 

0.92 

0.35 

0.80 

1 

850.6 

SC 

1.88s 

0.400s 

0.92 

0.35 

0.80 

1 

830.0 

ML 

1.457s 

0.464s 

0.94 

0.40 

0.35 

0.10 

2 

845. 7V 

CL 

1.589s 

0.422s 

0.97 

0.23 

0.05 

0.10 

2 

854.1 7 

OL 

1.165s 

0.533s 

0.955 

0.265 

0.20 

0.25 

2 

CH 

1.517  s 

0.457s 

0.98 

0.30 

0.05 

0.10 

2 

MH 

1.0601 

0.5471 

0.94 

0.30 

0.35 

0.10 

2 

830.0 

OH 

0.841 1 

0.8924 

0.955 

0.265 

0.20 

0.25 

2 

866. 77 

PT 

0.25 

0.70 

0.92 

0.40 

0.05 

0.50 

2 

830.0 

SMSC 

(MC) 

1 .60 1 

0.3961 

0.92 

0.40 

0.80 

0.05 

1 

830.0 

CLML 

(CM) 

1.617  s 

0.397s 

0.96 

0.30 

0.20 

0.10 

2 

830.0 

EV 

1.876s 

0.320s 

0.92 

0.40 

0.80 

0.00 

1 

830.0 

CO 

2.185 

0.020 

0.906 

0.406 

0.00 

850.0 

AS 

2.500 

0.020 

0.946 

0.1256 

0.00 

880.0 

RO 

2.700 

0.020 

0.89 

0.40 

0.00 

800.0 

SN 

0.920 

0.020 

0.90 

0.70 

0.00 

1  =  coarse,  2  =  fine 

1  Sullivan  et  al.  (1997),  p.  28,  38,  and  51 

2  Tarnawski  et  al.  (1997),  p.  96 

3  Farouki  (1981),  p.  136 

4  Yu  et  al.  (1993) 

5  Guymon,  G.L.  et  al.  (1993),  p.  21,  25 

6  Balick,  L.K.  et  al.  (1981),  Table  B1 

7  Peck,  L.  (2002),  Table  21 
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Table  9.2.2  FASST  Default  Soil 


Moisture  Properties. 


Soil 

Type 

Saturated 

Hydraulic 

Conductivity 

Ksat 

(cm/s) 

Residual 

Water 

Content 

0r 

(m'/m3) 

Maximum 

Water 

Content 

0max 

(m3/m3) 

Van 

Genuchten 

exponent 

nVG 

Bubbling 

Pressure 

Head 

CtvG 

(cm) 

Cl 

Coeff. 

Cl 1 

Cl 

Coeff. 

Cl 

RCI 

Coeff. 

Cl 

RCI 

Coeff. 

Cl 

GW 

0.000262 

0.01 

0.2962 

1.53 

22.6* 

GP 

0.00010562 

0.01 

0.4032 

2.23 

22.6* 

GM 

0.0000672 

0.01 

0.3242 

1.23 

32.7* 

GC 

0.000013891 

0.01 

0.344 

1.53 

23.2* 

SW 

0.000023612 

0.01 

0.3202 

1.253 

38.5473 

3.987 

0.81500 

3.987 

0.8150 

SP 

0.000007412 

0.01 

0.4152 

2.53 

38.5473 

3.987 

0.81500 

3.987 

0.8150 

SM 

7.987e-062 

0.01 

0.5262 

1.43 

68.5473 

8.749 

-1.1949 

12.542 

-2.955 

SC 

0.000000292 

0.01 

0.4002 

1.53 

58.5473 

9.056 

-1.3566 

12.542 

-2.955 

ML 

0.000000572 

0.01 

0.4642 

1.53 

33.9* 

10.225 

-1.565 

11.936 

-2.407 

CL 

7.7e-082 

0.01 

0.4222 

1.343 

53.5* 

10.998 

-1.848 

15.506 

-3.530 

OL 

0.000097222 

0.01 

0.5332 

1.343 

28.6* 

10.977 

-1.754 

17.399 

-3.584 

CH 

4.8e-082 

0.01 

0.4572 

1.53 

32.9* 

13.816 

-5.583 

13.686 

-2.705 

MH 

1.56-08* 

0.01 

0.5472 

1.343 

39.0* 

12.321 

-2.044 

23.641 

-5.191 

OH 

0.000883* 

0.01 

0.8924 

1.343 

29.3* 

13.046 

-2.172 

12.189 

-1.942 

PT 

0.000014 

0.15 

0.70 

1.34 

38.6 

SMSC 

(MC) 

7.7e-072 

0.01 

0.3962 

1.53 

23.5473 

9.056 

-1.3566 

12.542 

-2.955 

CLML 

(CM) 

1.26e-072 

0.01 

0.3972 

1.343 

32.9* 

9.454 

-1.3850 

14.236 

-3.137 

EV 

0.000023612 

0.01 

0.3202 

1.253 

38.5473 

3.987 

0.81500 

3.987 

0.8150 

CO, 

AS, 

RO, 

SN 

0.0 

0.001 

0.02 

1  Sullivan  et  al.  (1997),  p.  28 

2  Guymon,  et  al.  (1993),  p.  21,  25,  47-59 

3  Jordan  (2000) 

4  Yu  et  al.  (1993) 


Parameters  in  Table  9.2.2  are  used  to  calculate  the  moisture  content  and  resulting  strength 
of  the  soil.  They  are  needed  for  the  soil  moisture  and  soil  strength  subroutines  of  FASST. 
Unlike  other  models,  the  soil  strength  requires  moisture  content  by  weight  and  not  by 
volume.  The  moisture/strength  equations  are  a  result  of  curve  fitting.  They  are 

Cl  =  exp[c,  +  c2  ln(MC)]  (9.4) 

RCI  =  exp[Cj  +  c2  ln(  MC)] 

where  Cl  is  the  cone  index,  RCI  is  the  rating  cone  index,  and  MC  is  the  moisture  content 
by  weight  in  percent  (Sullivan  et  al.  1997).  The  rating  cone  index  is  used  to  quantify  the 
strength  of  soils  after  traffic  has  passed. 

The  relationship  between  volumetric  moisture  content  and  pressure  head  is  highly 
nonlinear.  Following  the  work  of  van  Genuchten  (1980)  it  is 
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(9.5) 


where  Q,  is  the  residual  volumetric  water  content,  6max  is  the  maximum  volumetric  water 
content,  aVG  (cm.)  is  a  constant  related  to  the  reciprocal  of  the  bubbling  pressure  head, 
nvc  is  a  constant  dependent  on  the  distribution  of  pores,  and  mvG  =  1  -  1  !nVQ.  The 
hydraulic  conductivity  also  depends  on  the  volumetric  moisture  content  such  that  (van 
Genuchten  1980) 


K  =  Ksalw' 


,0.5 


l-(l  -wUm'a)" 


n2 


(9.6) 


where 


w  = 


0-0r 

^max  ~  0r 


is  the  relative  soil  moisture. 


(9.7) 


9.3  Initial  Conditions 

FASST  requires  very  little  infonnation  concerning  the  initial  conditions.  At  a  minimum, 
it  needs  to  know  the  initial  snow  depth  and/or  ice  thickness  on  the  surface.  If  the  initial 
soil  temperature  and  moisture  profiles  are  known,  the  user  can  also  input  this 
information. 


9.4  Site  Information 

Several  parameters  describing  the  site  are  needed  by  FASST.  They  are 

Site  Latitude  (+  =  North) 

Site  Longitude  (+  =  East  from  Prime  Meridian) 

Site  Elevation  (m,ft) 

Time  Offset  between  Local  and  Universal  Time  Conversion  (UTC)  (hr) 

[local  -  UTC] 

Slope  (Degrees  from  horizontal) 

Aspect  Angle  (Degrees  from  North,  +  =  clockwise) 

Surface  Vegetation  0  =  unknown 

8  =  grass/pasture/steppe/meadow 
50  =  mixed  trees 
999  =  other 

The  vegetation  enumerations  are  based  on  SEDRIS  nomenclature.  At  present,  the 
vegetation  simply  changes  the  surface  albedo  and  emissivity.  The  next  release  of  the 
model  will  include  a  two-layer  vegetation  model  including  root  uptake  of  soil  moisture. 
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Chapter  10 

Model  Initialization 


10.0  Introduction 

Before  any  calculations  can  be  made,  the  surface  conditions  as  well  as  the  temperature, 
volumetric  moisture  content,  and  soil  properties  profiles  must  be  detennined.  The  initial 
nodal  spacing  and  number  of  nodes  are  also  calculated. 


10.1  Surface  Conditions 

The  presence  of  snow  or  ice  affects  both  the  soil  temperature  and  the  volumetric  moisture 
content.  It  also  influences  the  radiative  properties  of  the  surface,  specifically  the  albedo 
and  emissivity. 

The  albedo  is  a  measure  of  the  amount  of  solar  radiation  that  is  reflected  by  the  surface.  It 
is  the  ratio  of  the  reflected  to  the  total  incoming  solar  radiation.  The  greater  the  albedo, 
the  less  shortwave  radiation  is  absorbed  by  the  surface.  Value  of  the  albedo  for  a  snow- 
covered  soil  range  between  as  =  0.78  or  0.55  following  Pomeroy  et  al.  (1998),  depending 
on  the  age  of  the  snow  (see  Chapter  7.2.7).  For  ice-covered  soil  a;  =  0.70.  The  albedo 
values  for  the  USCS  soil  types  are  given  in  Table  9.2.1. 

The  emissivity  of  a  material  is  the  ratio  of  its  emitted  energy  to  emitted  infrared  radiation 
in  the  spectral  band  from  xxx  to  yyy  to  that  of  a  black  body,  or  perfect  absorber,  at  the 
same  temperature.  It  influences  the  amount  of  longwave  radiation  that  the  surface  loses 
to  the  atmosphere.  If  the  surface  is  snow-  or  ice-covered,  the  emissivity  is  ss  =  0.98  or 
a;  =  0.90,  respectively.  The  emissivity  values  for  the  USCS  soils  are  given  in  Table  9.2.1. 


10.2  Temperature  and  Volumetric  Moisture  and  Ice  Content  Profiles 

If  the  user  has  any  initial  soil  temperature  or  volumetric  moisture  content  data,  FASST 
incorporates  these  into  the  calculated  profiles.  For  both  parameters  FASST  assumes  a 
linear  profile  as  a  function  of  depth  between  any  consecutive  fixed  values. 

If  no  initial  surface  temperature  is  given,  FASST  assumes  that  the  soil  surface 
temperature  is  equal  to  the  air  temperature  if  there  is  no  snow  or  ice  on  the  ground.  If 
snow  or  ice  is  present,  FASST  sets  the  soil  surface  temperature  to  0.0°C.  If  the 
temperature  of  the  lowest  model  layer  is  not  provided,  FASST  assigns  a  value  equal  to 
the  surface  temperature  or  the  lowest  given  nodal  temperature. 

FASST  presumes  that  the  initial  volumetric  soil  moisture  content  at  the  surface  is  equal  to 
the  residual  moisture  allowed  plus  30%  of  the  difference  between  the  maximum  and 
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residual  values  unless  the  user  provides  a  value.  The  bottom  soil  moisture  is  assigned  the 
average  value  of  the  residual  and  maximum  moisture  contents  for  the  USCS  soil  type 
present  in  Table  9.2.2. 

If  the  temperature  of  the  soil  is  frozen  at  any  depth,  the  volumetric  ice  content  at  that 
position  equals  0.95*(soil  moisture)* (water  density)/(ice  density). 


10.3  Node  Position  and  Number 

FASST  requires  a  minimum  soil  depth  of  1.0  m.  If  the  user-specified  depth  is  less  than 
this,  FASST  extends  the  bottom  soil  layer  until  the  total  thickness  is  1.0  m. 

Node  spacing  is  detennined  by  both  the  number  of  soil  temperature  and  moisture 
observations  as  well  as  the  number  of  soil  layers.  If  measured  temperature  and/or 
moisture  values  exist,  a  node  is  placed  at  each  observation  location.  Nodes  also  occur  at 
any  layer  boundary  if  there  are  different  soil  types  as  well  as  the  soil  surface  and  bottom. 
If  no  initial  measurements  are  provided,  10  nodes  are  equally  spaced  through  the  total 
soil  thickness  and  at  any  layer  boundaries,  the  surface,  and  bottom.  Soil  type  dependent, 
temperature  and  moisture  independent,  parameters  are  assigned  to  the  nodes  depending 
on  which  layer  they  are  in.  Nodes  located  at  layer  boundaries  have  the  properties  of  the 
upper  layer. 


10.4  Time  Step  Determination 


Once  the  soil  type,  temperature,  and  volumetric  moisture  content  are  known  at  a  node, 
the  corresponding  hydraulic  and  thermal  properties  of  the  soil  are  calculated.  The  thermal 
model  calculation  time  step,  At  (sec),  is  then  calculated  based  on  the  numerical  stability 
requirement  for  the  soil  temperature.  Namely  for  each  node, 


At 


xth,i 


A,)3 


<0.5 


At  =  min 


0.45  (Az;) 


(10-1) 


2 

where  ktihi  (rn  /s)  is  the  thennal  diffusivity  of  the  soil  and  Azj(m)  is  the  nodal  spacing  and 
the  subscript  i  refers  to  the  node.  Since  the  thermal  diffusivity  varies  with  soil  type  and 
temperature,  the  time  step  also  varies.  As  indicated  in  Equation  (10.1),  the  minimum 
value  needed  to  satisfy  all  nodal  conditions  is  used. 

With  the  determination  of  the  number  of  nodes  and  their  placement,  as  well  as  the 
assignment  of  nodal  properties,  initialization  of  FASST  is  complete. 
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Chapter  11 

Model  Outputs 


11.0  Introduction 

FASST  produces  two  output  flies.  One  contains  information  pertaining  to  the  surface 
conditions  and  the  other  contains  soil  profile  infonnation.  Details  concerning  the  actual 
setup  of  the  two  files  may  be  found  in  the  FASST  User’s  Manual  (Frankenstein  2004). 


11.1  Surface  Condition  Information 

At  each  time  step  the  following  meteorological  data  and  surface  data  are  output. 

Year 

Day  of  the  Year 
Hour  of  the  Day 
Minute  of  the  Hour 
Air  Pressure  (mbar) 

Air  Temperature  ( K) 

Relative  Humidity  (%) 

Wind  Speed  (m/s) 

Wind  Direction  (degrees  from  N,  +  =  clockwise) 

Precipitation  Rate  (mm/hr) 

Precipitation  Type  (0  =  unknown,  1  =  none,  2  =  rain,  3  =  snow,  4  =  freezing  rain) 
Low  Cloud  Amount  ( tenths ) 

Low  Cloud  Height  (km) 

Low  Cloud  Type 

(If  BTRA  (VITD):0  =  none,  1  =  stratus,  2  =  stratocumulus,  3=cumulus 
Else  (SEDRIS):  0  =  none 

1=  cumulus  humulis  or  cumulus  fractus 

2  =  cumulus  mediocris  or  congestus 

3  =  cumulonimbus  calvus  with/out  cumulus 

stratocumulus  or  stratus 

4  =  stratocumulus  cumulogentius 

5  =  other  stratocumulus  types 

6  =  stratus  nebulosus  and/or  startus  fractus 

7  =  startus  fractus  and/or  cumulus  fractus  of  bad  weather 

8  =  cumulus  and  stratocumulus 

9  =  cumulonimbus  capillatus 
99  =  no  clouds  visible 

Middle  Cloud  Amount  (tenths) 

Middle  Cloud  Height  (km) 
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Middle  Cloud  Type 

(If  BTRA  (VITD):  0  =  none,  1  =  altostratus,  2  =  altocumulus, 

3  =  altocumulus  castellanus 
Else  (SEDRIS):  0  =  none 

1  =  altostratus  translucidus 

2  =  altostratus  opacus  or  nimbostratus 

3  =  altocumulus  translucidus,  1  level 

4  =  altocumulus  translucidus,  many  levels,  varying 

5  =  altocumulus  translucidus  in  bands 

6  =  altocumulus  cumulogentis  or  cumulonimbo  genitus 

7  =  altocumulus  translucidus  or  opacus,  multi  layer 

8  =  altocumulus  castellanus  or  floccus 

9  =  chaotic  altocumulus 

99  =  no  clouds  visible 

High  Cloud  Amount  ( tenths ) 

High  Cloud  Height  (km) 

High  Cloud  Type 

(If  BTRA  (VITD):0  =  none,  1  =  cirrus,  2  =  cirrostratus,  3  =  cirrocumulus 
Else  (SEDRIS):  0  =  none, 

1  =  cirrus  fibratus, 

2  =  cirrus  spissatu,  patchy 

3  =  cirrus  spissatus  cumulonimbo  genitus 

4  =  cirrus  unicinus  and/or  fibratus 

5  =  cirrus  and/or  cirrostratus  <  45  above  horizon 

6  =  cirrus  and/or  cirrostratus  >  45  above  horizon 

7  =  cirrostratus  full  cover 

8  =  cirrostratus  not  full  cover 

9  =  cirrocumulus 

99  =  no  clouds  visible 

Total  Incoming  Solar  Radiation  (Direct  +  Diffuse)  ( W/m  ) 

Direct  Solar  Radiation  (W/m2) 

Diffuse  Solar  Radiation  (W/m2) 

Reflected  Solar  Radiation  (W/m2) 

Net  Solar  Radiation  (Total  -Reflected)  (J/m2) 

Incoming  Infrared  Radiation  (W/m~) 

Emitted  Infrared  Radiation  ( W/m  ) 

Net  Infrared  Radiation  (J/m2) 

Solar  Zenith  Angle  (degrees) 

Solar  Azimuth  Angle  (degrees) 

Snow  Depth  (m) 

Ice  Thickness  (m) 

Soil  Surface  Temperature  (K) 

Soil  Surface  Moisture  (0  -  5cm)  (vol/vol) 

Surface  Cone  Index 
Surface  Rating  Cone  Index 
Freeze/Thaw  Depth  (m) 
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State  of  Freeze/Thaw  Depth  (‘f  rozen/’t’hawed/’u’nfrozen) 

Visibility  (km) 

Aerosol  Type  (0  =  unknown,  1  =  none,  2  =  mist,  3  =  dust,  4  =  smoke,  5  =  haze, 
6  =  ocean  spray,  7  =  sand,  8  =  volcanic  ash,  9  =  other) 

Snow  density  ( g/cm 3) 

Slippery  Indicator  for  mobility  (0  =  not  slippery,  1  =  wet,  2  =  ice,  3  =  snow) 


11.2  Soil  Profile  Information 

At  each  time  step,  the  following  information  is  output  for  each  node.  Node  1  is  at  the 
surface. 

Year 

Day  of  the  Year 
Hour  of  the  Day 
Minute  of  the  Hour 
Node  Number 
USCS  Soil  Type 
Soil  Temperature  ( K) 

Soil  Moisture  ( vol/vol ) 

Ice  Content  ( vol/vol) 

Total  Moisture  (Soil  Moisture  +  Ice  Content)  ( vol/vol) 

Frozen/Thawed  State 


11.3  Input/Output  Data  Source  Information 

Every  user  of  FASST  will  not  have  all  of  the  input  data  needed.  As  a  result,  FASST 
either  uses  default  values  or  calculates  the  parameters  internally.  If  a  default  value  is 
employed,  a  flag  of  “1”  is  written  to  this  file  along  with  the  variable  name  and  the 
corresponding  value.  If  the  model  calculates  the  parameter,  a  flag  of  “2”  is  written  to  the 
file  along  with  the  name  of  the  variable.  An  example  of  such  a  file  for  the  Grayling,  MI, 
input  data  is  shown  in  Figure  5  of  the  FASST  User’s  Manual  (Frankenstein  2004).  As  a 
default  this  file  is  named  “inferred. out.” 


11.4  References 
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